Table of contents

  • Importing Libraries
  • Configuring Visualization Parameters
  • Configuring Other Notebook Parameters
  • Pre-installing Custom Functions
  • Practicing in Stages
    • DICOM in Python
      • Extraction of DICOM Files
      • Visualization of 2D DICOM Files
      • Visualization of 3D Data Stored as Multiple 2D DICOM Files
      • Extraction of DICOM Files by Auto-Detection
    • NIfTI in Python
      • Extraction of the NIfTI File
      • Visualization of Extracted Images from the NIfTI File
      • Writing of the NIfTI file
    • Preprocessing in Python
      • Orientation
      • Canonicalization
      • Resampling

Importing Libraries¶

In [1]:
# Import Pydicom before using Pydicom functions

# DICOM (Digital Imaging and Communications in Medicine) is the standard protocol for
# the management and transmission of medical images and related data, and is used by
# many healthcare facilities

# Pydicom is a pure Python package for working with DICOM files
import pydicom
In [2]:
# The `dicom2nifti` library converts DICOM files to the NIfTI file and supports most anatomical
# CT and MRI data, especially MRI, which supports most 4D data such as DTI data and fMRI data
import dicom2nifti

# NiBabel is a pure Python package that supports a growing number of neuroimaging file formats,
# including the NIfTI1 format and the NIfTI2 format, which is an extension of the former

# NiBabel provides both format-independent access to high-level neuroimaging and format-specific APIs with
# different levels of access to all available information in a particular file format

# NiBabel's API provides full or selective access to file header information (metadata), while image data
# is provided via NumPy arrays
import nibabel as nib
import nibabel.processing
In [3]:
# The `pathlib` module is similar to the `os.path` module, but `pathlib` provides a more
# advanced and convenient interface than `os.path`

# It is possible to use `pathlib` to represent file paths as specialized `Path` objects
# instead of plain strings
from pathlib import Path

# SimpleITK is a simplified programming interface to the algorithms and data structures
# of the Insight Toolkit (ITK) for segmentation, registration and advanced image analysis

# SimpleITK supports bindings for multiple programming languages including C++, Python, R,
# Java, C#, Lua, Ruby and TCL

# Combining SimpleITK’s Python bindings with the Jupyter notebook web application creates
# an environment which facilitates collaborative development of biomedical image analysis
# workflows

# In this project, SimpleITK provides the ability to automatically detect and read all DICOM
# files so that the user does not have to manage file reading or slice sorting
import SimpleITK as sitk
In [4]:
import numpy as np

import torch
from torchvision.utils import make_grid
In [5]:
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
In [6]:
from datetime import datetime
from functools import wraps
import itertools
import math
import random
import reprlib
from scipy import ndimage
import sys

Configuring Visualization Parameters¶

In [7]:
%matplotlib inline
In [8]:
XINHUI = "#7a7374"
XUEBAI = "#fffef9"
YINBAI = "#f1f0ed"
YINHUI = "#918072"

figure_size = (16, 9)
In [9]:
custom_params = {
    "axes.axisbelow": True,
    "axes.edgecolor": YINBAI,
    "axes.facecolor": XUEBAI,
    "axes.grid": True,
    "axes.labelcolor": XINHUI,
    "axes.spines.right": False,
    "axes.spines.top": False,
    "axes.titlecolor": XINHUI,
    "figure.edgecolor": YINBAI,
    "figure.facecolor": XUEBAI,
    "grid.alpha": 0.8,
    "grid.color": YINBAI,
    "grid.linestyle": "--",
    "grid.linewidth": 1.2,
    "legend.edgecolor": YINHUI,
    "patch.edgecolor": XUEBAI,
    "patch.force_edgecolor": True,
    "text.color": XINHUI,
    "xtick.color": YINHUI,
    "ytick.color": YINHUI,
}

mpl.rcParams.update(custom_params)

Configuring Other Notebook Parameters¶

In [10]:
reprlib_rules = reprlib.Repr()
reprlib_rules.maxother = 250

Pre-installing Custom Functions¶

In [11]:
sys.path.append("../")
In [12]:
from Modules import *

Practicing in Stages¶

DICOM in Python¶

Extraction of DICOM Files¶

In [13]:
# The dataset used in this practice project is a very small subset of CT images extracted from
# the Cancer Imaging Archive (TCIA), which contains intermediate slices of all CT images
# in which valid age, mode, and contrast markers can be found

# This dataset is provided by the Kaggel Datasets called CT Medical Imaging
# (https://www.kaggle.com/datasets/kmader/siim-medical-images),
# license type is Attribution 3.0 Unported (CC BY 3.0)
# https://creativecommons.org/licenses/by/3.0/legalcode
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
# The main function of Pydicom to read and parse DICOM files is `read_file`
dicom_file = pydicom.read_file(dir_path + sample_dcm)

tabulation = Form_Generator()
tabulation.heading_printer("Initial understanding of DICOM file")

statements = [
    """
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
dicom_file = pydicom.read_file(dir_path + sample_dcm)
"""
]
tabulation.statement_generator(statements)

variables = ["dicom_file"]
values = [str(reprlib_rules.repr(dicom_file))]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "dicom_file[0x0028, 0x0010]",
    "dicom_file[0x0028, 0x0011]",
    "dicom_file[0x0018, 0x0015]",
    "dicom_file.Rows",
    "dicom_file.Columns",
    "dicom_file.BodyPartExamined",
    "dicom_file.keys()",
    "dicom_file.values()",
    "dicom_file.dir()",
    "dicom_file.dir('Image')",
]
results = [
    str(dicom_file[0x0028, 0x0010]),
    str(dicom_file[0x0028, 0x0011]),
    str(dicom_file[0x0018, 0x0015]),
    str(dicom_file.Rows),
    str(dicom_file.Columns),
    str(dicom_file.BodyPartExamined),
    str(reprlib_rules.repr(dicom_file.keys())),
    str(reprlib_rules.repr(dicom_file.values())),
    str(reprlib_rules.repr(dicom_file.dir())),
    str(reprlib_rules.repr(dicom_file.dir("Image"))),
]
tabulation.expression_generator(expressions, results, 1)
Initial understanding of DICOM file

    +-------------------------------------------------------+
    | Statement                                             |
    +-------------------------------------------------------+
    | dir_path = "../Datasets/Kaggle - CT Medical           |
    |     Images/dicom_dir/"                                |
    | sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"     |
    | dicom_file = pydicom.read_file(dir_path + sample_dcm) |
    +-------------------------------------------------------+
    +------------+------------------------------------------------+
    | Variable   | Value                                          |
    +------------+------------------------------------------------+
    | dicom_file | Dataset.file_meta                              |
    |            |  -------------------------------               |
    |            | (0002, 0000) File Meta Information Group       |
    |            |  Length  UL: 194                               |
    |            | (0002, 0001) Fil...Group Length                |
    |            |          UL: 524296                            |
    |            | (7fe0, 0010) Pixel Data                        |
    |            |    OW: Array of 524288 elements                |
    +------------+------------------------------------------------+
    +-----------------------------+-------------------------------+
    | Expression                  | Result                        |
    +-----------------------------+-------------------------------+
    | dicom_file[0x0028, 0x0010]  | (0028, 0010) Rows             |
    |                             |                     US: 512   |
    | dicom_file[0x0028, 0x0011]  | (0028, 0011) Columns          |
    |                             |                     US: 512   |
    | dicom_file[0x0018, 0x0015]  | (0018, 0015) Body Part        |
    |                             |  Examined                     |
    |                             |  CS: 'CHEST'                  |
    | dicom_file.Rows             | 512                           |
    | dicom_file.Columns          | 512                           |
    | dicom_file.BodyPartExamined | CHEST                         |
    | dicom_file.keys()           | dict_keys([(0008, 0000),      |
    |                             |  (0008, 0005), (0008, 0008),  |
    |                             |  (0008, 0016), (0008, 0018),  |
    |                             |  (0008, 0020), (0008, 0021),  |
    |                             |  (0008, 0022), ...095, 0010), |
    |                             |  (0097, 0000), (0097, 0010),  |
    |                             |  (0099, 0000), (0099, 0010),  |
    |                             |  (7003, 0000), (7003, 0010),  |
    |                             |  (7fe0, 0000), (7fe0, 0010)]) |
    | dicom_file.values()         | dict_values([(0008, 0000)     |
    |                             |  Group Length                 |
    |                             |         UL: 430, (0008, 0005) |
    |                             |  Specific Character Set       |
    |                             |         CS:...up Length       |
    |                             |                   UL: 524296, |
    |                             |  (7fe0, 0010) Pixel Data      |
    |                             |                      OW:      |
    |                             |  Array of 524288 elements])   |
    | dicom_file.dir()            | ['AccessionNumber',           |
    |                             |  'AcquisitionDate',           |
    |                             |  'AcquisitionNumber',         |
    |                             |  'AcquisitionTime',           |
    |                             |  'BitsAllocated',             |
    |                             |  'BitsStored', ...]           |
    | dicom_file.dir('Image')     | ['ImageDimensions',           |
    |                             |  'ImageFormat',               |
    |                             |  'ImageGeometryType',         |
    |                             |  'ImageLocation',             |
    |                             |  'ImageOrientation',          |
    |                             |  'ImageOrientationPatient',   |
    |                             |  ...]                         |
    +-----------------------------+-------------------------------+
In [14]:
# DICOM data has an attribute called `pixel_array` that provides more useful pixel data
# for uncompressed images that can be passed to the graphics library for viewing

# To use this attribute, the system must have the NumPy numeric package installed,
# since `pixel_array` returns a NumPy array
ct = dicom_file.pixel_array

tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from DICOM file")

statements = ["ct = dicom_file.pixel_array"]
tabulation.statement_generator(statements)

variables = ["ct"]
values = [str(reprlib_rules.repr(ct))]
tabulation.variable_generator(variables, values)

expressions = ["ct.shape"]
results = [str(ct.shape)]
tabulation.expression_generator(expressions, results)
Getting pixel data from DICOM file

    +-----------------------------+
    | Statement                   |
    +-----------------------------+
    | ct = dicom_file.pixel_array |
    +-----------------------------+
    +----------+------------------------------------------------+
    | Variable | Value                                          |
    +----------+------------------------------------------------+
    | ct       | array([[0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        ...,                                    |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0],                |
    |          |        [0, 0, 0, ..., 0, 0, 0]], dtype=uint16) |
    +----------+------------------------------------------------+
    +------------+------------+
    | Expression | Result     |
    +------------+------------+
    | ct.shape   | (512, 512) |
    +------------+------------+

Visualization of 2D DICOM Files¶

In [15]:
def image_display(image, ax, title, cmap):
    ax.imshow(image, cmap)
    ax.grid(False)
    ax.set_title(title, loc="center", pad=10)
    x_ticks = list(range(0, image.shape[1], 50))
    y_ticks = list(range(0, image.shape[0], 50))
    ax.set(xticks=x_ticks, xticklabels=x_ticks,
           yticks=y_ticks, yticklabels=y_ticks)
    ax.set_xlim(left=0)
    ax.set_ylim(top=0)
    ax.minorticks_on()
    return ax


# Path classes are divided between pure paths, which provide purely computational
# operations without I/O, and concrete paths, which inherit from pure paths but also
# provide I/O operations
path_object = Path(dir_path)
# `PurePath.name` returns a string representing the final path component, excluding
# the drive and root directory (if any)

# When a path points to a directory, `Path.iterdir` generates a path object of the directory's
# contents

# `Path.is_file` returns True if the path points to a normal file or to a symbolic link
# to a normal file, False if it points to another type of file
random_dicom_path = random.choice(
    [
        file.name
        for file in path_object.iterdir()
        if file.is_file() & (pydicom.read_file(file).BodyPartExamined != "CHEST")
    ]
)

random_dicom_file = pydicom.read_file(dir_path + random_dicom_path)

plt.rcParams["figure.figsize"] = (
    figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)

fig, axs = plt.subplots(nrows=2, ncols=2)

image_display(
    ct,
    axs[0, 0],
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
    cmap="gray",
)

image_display(
    ct,
    axs[0, 1],
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
    cmap="bone",
)

image_display(
    random_dicom_file.pixel_array,
    axs[1, 0],
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
    cmap="gray",
)

image_display(
    random_dicom_file.pixel_array,
    axs[1, 1],
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
    cmap="bone",
)


fig.suptitle(
    "Visual Comparison of CT Medical Image Display Using Different Colormaps",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [16]:
def annotation_color(cmap):
    cmap_color = mpl.colormaps[cmap]
    return cmap_color(0.85)


def first_max(list):
    max_element = list[0]
    for element in list[1:]:
        if len(element) > len(max_element):
            max_element = element
    return max_element


def text_aligner(text):
    line_list = text.split("\n")
    max_length = len(first_max(line_list))
    line_list = [
        (":" + " " * (max_length - len(line))).join(line.split(":"))
        if len(line) < max_length
        else line
        for line in line_list
    ]
    text = "\n".join(line_list)
    return text


def plane_information_annotator(func):
    @wraps(func)
    def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
        plane = func(ax, cmap, dicom_file, fontsize_adjustment)
        text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}"
        text += f"\nImage Plane: {plane}"
        text += f"\nPatient Position: {dicom_file.PatientPosition}"
        text += f"\nSlice Thickness: {dicom_file.SliceThickness:.1f} mm"
        text += f"\nSlice Location: {dicom_file.SliceLocation:.1f} mm"
        ax.text(
            0.575,
            0.825,
            text_aligner(text),
            transform=ax.transAxes,
            fontsize=9 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )

    return wrapper


def orientation_annotator(func):
    @wraps(func)
    def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
        plane, top, right, bottom, left = func(dicom_file)
        ax.text(
            0.4875,
            0.95,
            top,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.4875,
            0.025,
            bottom,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.95,
            0.4875,
            right,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        ax.text(
            0.025,
            0.4875,
            left,
            transform=ax.transAxes,
            fontweight="heavy",
            fontsize=11 + fontsize_adjustment,
            fontfamily="monospace",
            c=annotation_color(cmap),
        )
        return plane

    return wrapper


def image_orientation(func):
    @wraps(func)
    def wrapper(*args, **kwargs):
        Left, Right, Anterior, Posterior, Superior, Inferior = (
            "L",
            "R",
            "A",
            "P",
            "S",
            "I",
        )
        axial_plane = [Anterior, Right, Posterior, Left]
        coronal_plane = [Superior, Left, Inferior, Right]
        sagittal_plane = [Superior, Anterior, Inferior, Posterior]
        plane = func(*args, **kwargs)
        # Patient Position (0018,5100) specifies the position of the patient relative to
        # the space of the imaging device and is used for annotation purposes only,
        # but does not provide a precise mathematical relationship between the patient and
        # the imaging device
        patient_position = dicom_file.PatientPosition
        if plane == "Axial":
            term_list = ["S", "DR", "P", "DL"]
            top = 4 - term_list.index(patient_position[2:])
            if patient_position.startswith("HF"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                axial_plane[index] if index < 4 else axial_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        elif plane == "Coronal":
            term_list = ["HF", "RF", "FF", "LF"]
            top = 4 - term_list.index(patient_position[:2])
            if patient_position.endswith("S"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                coronal_plane[index] if index < 4 else coronal_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        elif plane == "Sagittal":
            term_list = ["HF", "PF", "FF", "AF"]
            top = 4 - term_list.index(patient_position[:2])
            if patient_position.endswith("DL"):
                right = top + 1
                left = right + 2
            else:
                left = top + 1
                right = left + 2
            bottom = top + 2
            top, right, bottom, left = [
                sagittal_plane[index] if index < 4 else sagittal_plane[index - 4]
                for index in [top, right, bottom, left]
            ]
        else:
            top, right, bottom, left = np.repeat("", 4)
        return plane, top, right, bottom, left

    return wrapper


@plane_information_annotator
@orientation_annotator
@image_orientation
def get_plane_information(dicom_file):
    # The DICOM standard defines a patient-oriented reference coordinate system (RCS) that
    # enables the user to measure the position and orientation of the image relative to
    # the patient

    # Image Orientation (Patient) (0020,0037) specifies the cosine of the orientation of
    # the first row and column relative to the patient, which should be provided as a pair,
    # where the row values for the x, y, and z axes are followed by the column values for
    # the x, y, and z axes

    # The orientation of the axes is determined solely by the orientation of the patient, i.e.,
    # the RCS and Image Orientation (Patient) (0020,0037) values specify the orientation of
    # the image frame rows and columns
    IOP = dicom_file.ImageOrientationPatient
    row_xyz = IOP[:3]
    column_xyz = IOP[3:]
    # `numpy.cross` returns the cross product of two vector arrays

    # The geometric interpretation of the cross product is that two vectors determine a plane,
    # and the cross product points in a direction different from both vectors

    # The Patient-Based Coordinate System is a right-handed system, i.e., the vector
    # cross product of a unit vector along the positive x-axis and a unit vector along
    # the positive y-axis is equal to a unit vector along the positive z-axis
    plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
    # The Image Type (0008,0008) identifies important image identification features
    # and it returns the Pixel Data Characteristics, Patient Examination Characteristics,
    # Modality Specific Characteristics, and Implementation Specific Identifiers

    # While the third value (Modality Specific Characteristics) should identify any
    # specialization specific to the Image Information Object Definition (IOD), it
    # and the values that follow it are not mandatory, unlike the first two values,
    # and therefore some DICOM data does not have this value
    if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
        return "Sagittal"
    elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
        return "Coronal"
    elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
        return "Axial"
    else:
        "Unknown"


def get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age):
    text = f"Patient ID: {dicom_file.PatientID}"
    text += f"\nPatient's Sex: {dicom_file.PatientSex}"
    try:
        patien_age = dicom_file.PatientAge.strip(str(0))[:-1]
    except AttributeError:
        patien_age = patien_age
    text += f"\nPatient's Age: {patien_age} y"
    text += f"\nModality: {dicom_file.Modality}"
    text += f"\nBody Part Examined: {dicom_file.BodyPartExamined}"
    ax.text(
        0.04,
        0.825,
        text_aligner(text),
        transform=ax.transAxes,
        fontsize=9 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def get_date_information(ax, cmap, dicom_file, fontsize_adjustment):
    date = datetime.strptime(dicom_file.StudyDate,
                             "%Y%m%d").strftime("%Y.%m.%d")
    text = f"Study Date: {date}"
    ax.text(
        0.04,
        0.05,
        text,
        transform=ax.transAxes,
        fontsize=7 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment):
    text = f"Manufacturer: {dicom_file.Manufacturer}"
    ax.text(
        0.525,
        0.05,
        text.rjust(36),
        transform=ax.transAxes,
        fontsize=7 + fontsize_adjustment,
        fontfamily="monospace",
        c=annotation_color(cmap),
    )


def DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment=0, patien_age="__"):
    ax.imshow(dicom_file.pixel_array, cmap)
    ax.grid(False)
    ax.set_title(title, loc="center", pad=10,
                 fontsize=12 + fontsize_adjustment)
    ax.set(xticks=[], yticks=[], frame_on=False)
    get_patient_information(ax, cmap, dicom_file,
                            fontsize_adjustment, patien_age)
    get_plane_information(ax, cmap, dicom_file, fontsize_adjustment)
    get_date_information(ax, cmap, dicom_file, fontsize_adjustment)
    get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment)


fig, axs = plt.subplots(2, 2, figsize=(
    figure_size[0] / 3 * 2, figure_size[1] / 4 * 5))

DICOM_2D_image(
    axs[0, 0],
    "gray",
    dicom_file,
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
)

DICOM_2D_image(
    axs[0, 1],
    "bone",
    dicom_file,
    f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
)

DICOM_2D_image(
    axs[1, 0],
    "gray",
    random_dicom_file,
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'gray'",
)

DICOM_2D_image(
    axs[1, 1],
    "bone",
    random_dicom_file,
    f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
    "colormap 'bone'",
)


fig.suptitle(
    "Visualization of CT Medical Image Display with Supplementary Information",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()

Visualization of 3D Data Stored as Multiple 2D DICOM Files¶

In [17]:
# The dataset used in this practice project is the MRI DICOM dataset of the head of
# a 52-year-old normal male mathematics professor

# The subject of the experiment is the author, who suffers from small vertical strabismus
# (hypertropia), which is a misalignment of the eyes, which can be seen in this dataset

# This dataset was provided by Lionheart, William R.B. in 2015, available online
# (https://zenodo.org/record/16956 or http://doi.org/10.5281/zenodo.16956),
# license type is Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)
# https://creativecommons.org/licenses/by-sa/4.0/legalcode
path_to_head_mri = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
# `Path.glob` selects the given relative pattern in the directory represented by the given path
# and yields all matching files (of any type)

# Since `Path.glob` returns a generator, it is converted to a list here for easy display
all_files = list(path_to_head_mri.glob("*"))

mri_data = []

for path in all_files:
    # Read the individual DICOM files one by one in the list returned by the path generator
    data = pydicom.read_file(path)
    mri_data.append(data)

tabulation = Form_Generator()
tabulation.heading_printer("Reading all DICOM files that match the specified pattern")

statements = [
    """
path_to_head_mri = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
all_files = list(path_to_head_mri.glob("*"))

mri_data = []

for path in all_files:
    data = pydicom.read_file(path)
    mri_data.append(data)
"""
]
tabulation.statement_generator(statements)

variables = ["path_to_head_mri", "all_files"]
values = [str(path_to_head_mri), str(reprlib_rules.repr(all_files))]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "all_files[0].parts",
    "all_files[0].parent",
    "all_files[0].name",
    "all_files[0].suffixes",
    "all_files[0].stem",
    "len(all_files)",
    "mri_data[0]",
    "len(mri_data)",
]
results = [
    str(all_files[0].parts),
    str(all_files[0].parent),
    str(all_files[0].name),
    str(all_files[0].suffixes),
    str(all_files[0].stem),
    str(len(all_files)),
    str(reprlib_rules.repr(mri_data[0])),
    str(len(mri_data)),
]
tabulation.expression_generator(expressions, results, 1)
Reading all DICOM files that match the specified pattern

    +---------------------------------------------------------+
    | Statement                                               |
    +---------------------------------------------------------+
    | path_to_head_mri = Path(                                |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a |
    |     Normal Male Human Aged 52/"                         |
    |     "ST000000/SE000001/"                                |
    | )                                                       |
    | all_files = list(path_to_head_mri.glob("*"))            |
    |                                                         |
    | mri_data = []                                           |
    |                                                         |
    | for path in all_files:                                  |
    |     data = pydicom.read_file(path)                      |
    |     mri_data.append(data)                               |
    +---------------------------------------------------------+
    +------------------+------------------------------------------+
    | Variable         | Value                                    |
    +------------------+------------------------------------------+
    | path_to_head_mri | ../Datasets/An MRI Dicom Data Set of the |
    |                  |  Head of a Normal Male Human Aged        |
    |                  |  52/ST000000/SE000001                    |
    | all_files        | [PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000000'),        |
    |                  |  PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000007'),        |
    |                  |  PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000009'),        |
    |                  |  PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000008'),        |
    |                  |  PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000006'),        |
    |                  |  PosixPath('../Datasets/An MRI Dicom     |
    |                  |  Data Set of the Head of a Normal Male   |
    |                  |  Human Aged                              |
    |                  |  52/ST000000/SE000001/MR000001'), ...]   |
    +------------------+------------------------------------------+
    +-----------------------+-------------------------------------+
    | Expression            | Result                              |
    +-----------------------+-------------------------------------+
    | all_files[0].parts    | ('..', 'Datasets', 'An MRI Dicom    |
    |                       |  Data Set of the Head of a Normal   |
    |                       |  Male Human Aged 52', 'ST000000',   |
    |                       |  'SE000001', 'MR000000')            |
    | all_files[0].parent   | ../Datasets/An MRI Dicom Data Set   |
    |                       |  of the Head of a Normal Male Human |
    |                       |  Aged 52/ST000000/SE000001          |
    | all_files[0].name     | MR000000                            |
    | all_files[0].suffixes | []                                  |
    | all_files[0].stem     | MR000000                            |
    | len(all_files)        | 27                                  |
    | mri_data[0]           | Dataset.file_meta                   |
    |                       |  -------------------------------    |
    |                       | (0002, 0000) File Meta Information  |
    |                       |  Group Length  UL: 214              |
    |                       | (0002, 0001) Fil...Group Length     |
    |                       |                     UL: 131084      |
    |                       | (7fe0, 0010) Pixel Data             |
    |                       |               OW: Array of 131072   |
    |                       |  elements                           |
    | len(mri_data)         | 27                                  |
    +-----------------------+-------------------------------------+
In [18]:
# DICOM files may not be sorted according to their actual image positions, this can
# be verified by checking Slice Location (0020,1041)

# Slice Location (0020,1041) identifies the relative position of the image plane in millimeters
unordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}

# For better viewing and processing of 3D data stored with multiple 2D DICOM files,
# these files must be sorted; otherwise the complete scanned image will become cluttered
# and useless
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)

ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}


tabulation = Form_Generator()
tabulation.heading_printer("Sorting of read DICOM files")

statements = [
    """
unordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}

mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)

ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
"""
]
tabulation.statement_generator(statements)

variables = ["unordered_slices", "ordered_slices"]
values = [
    str(unordered_slices),
    str(ordered_slices),
]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "mri_data_ordered[0]",
    "len(mri_data_ordered)",
    "len(unordered_slices)",
    "len(ordered_slices)",
]
results = [
    str(reprlib_rules.repr(mri_data_ordered[0])),
    str(len(mri_data_ordered)),
    str(len(unordered_slices)),
    str(len(ordered_slices)),
]
tabulation.expression_generator(expressions, results, 1)
Sorting of read DICOM files

    +-----------------------------------------------------------+
    | Statement                                                 |
    +-----------------------------------------------------------+
    | unordered_slices = {                                      |
    |     index: float(slice.SliceLocation) for index, slice in |
    |     enumerate(mri_data)                                   |
    | }                                                         |
    |                                                           |
    | mri_data_ordered = sorted(mri_data, key=lambda slice:     |
    |     slice.SliceLocation)                                  |
    |                                                           |
    | ordered_slices = {                                        |
    |     index: float(slice.SliceLocation) for index, slice in |
    |     enumerate(mri_data_ordered)                           |
    | }                                                         |
    +-----------------------------------------------------------+
    +------------------+------------------------------------------+
    | Variable         | Value                                    |
    +------------------+------------------------------------------+
    | unordered_slices | {0: 0.0, 1: 41.9999963629367, 2:         |
    |                  |  53.9999958207213, 3: 47.9999970362677,  |
    |                  |  4: 35.9999959546749, 5:                 |
    |                  |  5.99999663091323, 6: 137.999998321624,  |
    |                  |  7: 143.999998928727, 8:                 |
    |                  |  71.9999961590453, 9: 89.9999955528687,  |
    |                  |  10: 83.9999967682912, 11:               |
    |                  |  77.999996227574, 12: 149.999999502083,  |
    |                  |  13: 131.999997780749, 14:               |
    |                  |  23.9999946081714, 15: 17.9999979772582, |
    |                  |  16: 11.9999973042441, 17:               |
    |                  |  29.9999952815023, 18: 107.999995419197, |
    |                  |  19: 119.999996566542, 20:               |
    |                  |  95.9999960937442, 21: 65.9999961939969, |
    |                  |  22: 59.9999962290673, 23:               |
    |                  |  101.999994745866, 24: 125.999997173645, |
    |                  |  25: 155.999992554172, 26:               |
    |                  |  113.999995959439}                       |
    | ordered_slices   | {0: 0.0, 1: 5.99999663091323, 2:         |
    |                  |  11.9999973042441, 3: 17.9999979772582,  |
    |                  |  4: 23.9999946081714, 5:                 |
    |                  |  29.9999952815023, 6: 35.9999959546749,  |
    |                  |  7: 41.9999963629367, 8:                 |
    |                  |  47.9999970362677, 9: 53.9999958207213,  |
    |                  |  10: 59.9999962290673, 11:               |
    |                  |  65.9999961939969, 12: 71.9999961590453, |
    |                  |  13: 77.999996227574, 14:                |
    |                  |  83.9999967682912, 15: 89.9999955528687, |
    |                  |  16: 95.9999960937442, 17:               |
    |                  |  101.999994745866, 18: 107.999995419197, |
    |                  |  19: 113.999995959439, 20:               |
    |                  |  119.999996566542, 21: 125.999997173645, |
    |                  |  22: 131.999997780749, 23:               |
    |                  |  137.999998321624, 24: 143.999998928727, |
    |                  |  25: 149.999999502083, 26:               |
    |                  |  155.999992554172}                       |
    +------------------+------------------------------------------+
    +-----------------------+-------------------------------------+
    | Expression            | Result                              |
    +-----------------------+-------------------------------------+
    | mri_data_ordered[0]   | Dataset.file_meta                   |
    |                       |  -------------------------------    |
    |                       | (0002, 0000) File Meta Information  |
    |                       |  Group Length  UL: 214              |
    |                       | (0002, 0001) Fil...Group Length     |
    |                       |                     UL: 131084      |
    |                       | (7fe0, 0010) Pixel Data             |
    |                       |               OW: Array of 131072   |
    |                       |  elements                           |
    | len(mri_data_ordered) | 27                                  |
    | len(unordered_slices) | 27                                  |
    | len(ordered_slices)   | 27                                  |
    +-----------------------+-------------------------------------+
In [19]:
# Fill the 3D array in a slice-per-slice manner
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)

tabulation = Form_Generator()
tabulation.heading_printer(
    "One-time extraction of the overall 3D pixel array for all slices"
)

statements = [
    """
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
"""
]
tabulation.statement_generator(statements)

variables = ["full_volume"]
values = [str(reprlib_rules.repr(full_volume))]
tabulation.variable_generator(variables, values)

expressions = ["full_volume.shape", "full_volume.dtype"]
results = [str(full_volume.shape), str(full_volume.dtype)]
tabulation.expression_generator(expressions, results)
One-time extraction of the overall 3D pixel array for all slices

    +-----------------------------------------------+
    | Statement                                     |
    +-----------------------------------------------+
    | full_volume = [slice.pixel_array for slice in |
    |     mri_data_ordered]                         |
    | full_volume = np.array(full_volume)           |
    +-----------------------------------------------+
    +-------------+------------------------------------+
    | Variable    | Value                              |
    +-------------+------------------------------------+
    | full_volume | array([[[0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         ...,                       |
    |             |         [0,...     ...,            |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0],   |
    |             |         [0, 0, 0, ..., 0, 0, 0]]], |
    |             |         dtype=uint16)              |
    +-------------+------------------------------------+
    +-------------------+----------------+
    | Expression        | Result         |
    +-------------------+----------------+
    | full_volume.shape | (27, 256, 256) |
    | full_volume.dtype | uint16         |
    +-------------------+----------------+
In [20]:
def sort_checker(list_1):
    list_2 = list_1[1:]
    if all(a <= b for a, b in zip(list_1, list_2)):
        return "ascending"
    elif all(a >= b for a, b in zip(list_1, list_2)):
        return "descending"
    elif all(a == b for a, b in zip(list_1, list_2)):
        return "identical"
    else:
        return "unordered"


def title_adder(func):
    @wraps(func)
    def wrapper(inputs, ax, dict_slices, title=None, **kwargs):
        is_sorted = sort_checker(list(dict_slices.values()))
        if title is None:
            title = f"DICOM image grid with {is_sorted} slice locations"
        func(inputs, ax, dict_slices, **kwargs)
        ax.set_title(
            title,
            loc="center",
            pad=10,
        )

    return wrapper


def grid_annotator(func):
    @wraps(func)
    def wrapper(inputs, ax, dict_slices, fontsize_adjustment=0, **kwargs):
        x_positions_1, x_positions_2, y_positions, cmap = func(inputs, ax, **kwargs)
        for (y, x_1), key in zip(
            itertools.product(y_positions, x_positions_1), dict_slices.keys()
        ):
            text = f"No.{(key + 1):>03}"
            ax.text(
                x_1,
                y,
                text,
                transform=ax.transData,
                fontsize=7 + fontsize_adjustment,
                fontfamily="monospace",
                c=annotation_color(cmap),
            )
        for (y, x_2), value in zip(
            itertools.product(y_positions, x_positions_2), dict_slices.values()
        ):
            text = f"{value:>5.1f} mm"
            ax.text(
                x_2,
                y,
                text,
                transform=ax.transData,
                fontsize=7 + fontsize_adjustment,
                fontfamily="monospace",
                c=annotation_color(cmap),
            )

    return wrapper


@title_adder
@grid_annotator
def grid_DICOM_2D_image(inputs, ax, n=None, row_size=9, pad_value=1):
    if n is None:
        n = inputs.shape[0]
    nrows = math.ceil(n / row_size)
    cmap = "bone"

    # The only types supported by `torch.Tensor` are float64, float32, float16, complex64,
    # complex128, int64, int32, int16, int8, uint8 and bool, not uint16
    inputs = torch.Tensor(inputs.astype(int)).unsqueeze(1)

    x_positions_1 = [
        math.ceil(
            x * inputs[:row_size].shape[3]
            + 0.075 * inputs[:row_size].shape[3]
            + pad_value * x
        )
        for x in range(inputs[:row_size].shape[0])
    ]
    x_positions_2 = [
        math.ceil(
            x * inputs[:row_size].shape[3]
            + 0.55 * inputs[:row_size].shape[3]
            + pad_value * x
        )
        for x in range(inputs[:row_size].shape[0])
    ]
    y_positions = [
        math.ceil(
            y * inputs[:row_size].shape[2]
            + 0.125 * inputs[:row_size].shape[2]
            + pad_value * y
        )
        for y in range(nrows)
    ]

    for i in range(nrows):
        if len(inputs) > row_size:
            row_images = inputs[:row_size]
            inputs = inputs[row_size:]
        else:
            row_images = inputs[:]
        # Unlike the images in the MNIST dataset, the DICOM file images do not take values
        # in the range (0, 1), so the parameter `normalize` needs to be set to True in order to
        # move their values proportionally to the range (0, 1)

        # By querying the source code of `vision/torchvision/utils.py`, it is clear that any
        # single-channel image is converted to a three-channel image, and all are represented
        # as tensors, as shown below:
        # > if tensor.dim() == 2:  # single image H x W
        # >     tensor = tensor.unsqueeze(0)
        # > if tensor.dim() == 3:  # single image
        # >     if tensor.size(0) == 1:  # if single-channel, convert to 3-channel
        # >         tensor = torch.cat((tensor, tensor, tensor), 0)
        # >     tensor = tensor.unsqueeze(0)
        # >
        # > if tensor.dim() == 4 and tensor.size(1) == 1:  # single-channel images
        # >     tensor = torch.cat((tensor, tensor, tensor), 1)
        grid_row = make_grid(
            row_images, nrow=row_size, normalize=True, pad_value=pad_value
        )
        scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
        if i == 0:
            if nrows != 1:
                scalar_image = scalar_row
            else:
                scalar_image = np.full(
                    (
                        scalar_row.shape[0],
                        math.ceil(scalar_row.shape[1] / inputs.shape[0] * row_size),
                    ),
                    np.nan,
                )
                scalar_image[:, : scalar_row.shape[1]] = scalar_row
        elif scalar_image.shape[1] == scalar_row.shape[1]:
            scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
        else:
            image_to_concat = np.full(
                (scalar_image.shape[0] + scalar_row.shape[0], scalar_image.shape[1]),
                np.nan,
            )
            image_to_concat[: scalar_image.shape[0], :] = scalar_image
            image_to_concat[scalar_image.shape[0] :, : scalar_row.shape[1]] = scalar_row
            scalar_image = image_to_concat
    ax.imshow(scalar_image, cmap=cmap)
    ax.set(xticks=[], yticks=[], frame_on=False)
    ax.grid(False)
    return x_positions_1, x_positions_2, y_positions, cmap


unordered_full_volume = [slice.pixel_array for slice in mri_data]
unordered_full_volume = np.array(unordered_full_volume)
descending_slices = {
    index: float(slice.SliceLocation)
    for index, slice in enumerate(mri_data_ordered[::-1])
}

fig, axs = plt.subplots(3, 1, figsize=(figure_size[0], figure_size[1] / 7 * 9))
grid_DICOM_2D_image(unordered_full_volume, axs[0], unordered_slices)
grid_DICOM_2D_image(full_volume, axs[1], ordered_slices)
grid_DICOM_2D_image(full_volume[::-1], axs[2], descending_slices)

fig.suptitle(
    "Visual Comparison of DICOM Image Grids with Unordered and Ordered Slice Locations",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [21]:
fig, axs = plt.subplots(9, 3, figsize=(
    figure_size[0] / 3 * 2, figure_size[1] / 2 * 7))

for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))):
    DICOM_2D_image(
        axs[i][j],
        "bone",
        mri_data_ordered[slice_counter],
        f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm "
        "slice location",
        fontsize_adjustment=-3,
        patien_age="52",
    )

fig.suptitle(
    "Visualization of MRI Medical Images with Information Displayed in Order of Slice Location",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()

Extraction of DICOM Files by Auto-Detection¶

In [22]:
# The path object must be converted to a string in order for the SimpleITK package to process
# the path information
path_string = str(path_to_head_mri)
# For some image formats such as DICOM, the images also contain associated metadata,
# which is not loaded by the reader by default for time-saving reasons

# `ImageSeriesReader` class provides the ability to read a series of image files into a
# SimpleITK image, and once a series of images has been read, the metadata can be accessed
# directly from the reader

# `GetGDCMSeriesIDs` provides the ability to get all seriesIDs from a DICOM dataset
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
# `GetGDCMSeriesFileNames` generates a series of filenames from the directory containing the
# DICOM dataset and series ID

# Note that the resulting series of filenames (which are actually file paths) are automatically
# sorted
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    path_string, series_ids[0]
)

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)

# Finally, the reader is executed by calling `Execute` to get the required data
image_data = series_reader.Execute()

tabulation = Form_Generator()
tabulation.heading_printer(
    "Automatic recognition and recall of series filenames for DICOM images"
)

statements = [
    """
path_string = str(path_to_head_mri)
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    path_string, series_ids[0]
)

series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)

image_data = series_reader.Execute()
"""
]
tabulation.statement_generator(statements)

variables = ["series_ids", "series_file_names", "series_reader", "image_data"]
values = [
    str(reprlib_rules.repr(series_ids)),
    str(reprlib_rules.repr(series_file_names)),
    str(reprlib_rules.repr(series_reader)),
    str(reprlib_rules.repr(image_data)),
]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "type(series_ids)",
    "len(series_ids)",
    "type(series_file_names)",
    "len(series_file_names)",
    "len(image_data)",
    "image_data.GetSize()",
]
results = [
    str(type(series_ids)),
    str(len(series_ids)),
    str(type(series_file_names)),
    str(len(series_file_names)),
    str(len(image_data)),
    str(image_data.GetSize()),
]
tabulation.expression_generator(expressions, results)
Automatic recognition and recall of series filenames for DICOM images

    +----------------------------------------------------------+
    | Statement                                                |
    +----------------------------------------------------------+
    | path_string = str(path_to_head_mri)                      |
    | series_ids =                                             |
    |     sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string) |
    | series_file_names =                                      |
    |     sitk.ImageSeriesReader.GetGDCMSeriesFileNames(       |
    |     path_string, series_ids[0]                           |
    | )                                                        |
    |                                                          |
    | series_reader = sitk.ImageSeriesReader()                 |
    | series_reader.SetFileNames(series_file_names)            |
    |                                                          |
    | image_data = series_reader.Execute()                     |
    +----------------------------------------------------------+
    +-------------------+-----------------------------------------+
    | Variable          | Value                                   |
    +-------------------+-----------------------------------------+
    | series_ids        | ('1.3.46.67058...1413262801702',)       |
    | series_file_names | ('../Datasets/...0001/MR000000',        |
    |                   |  '../Datasets/...0001/MR000001',        |
    |                   |  '../Datasets/...0001/MR000002',        |
    |                   |  '../Datasets/...0001/MR000003',        |
    |                   |  '../Datasets/...0001/MR000004',        |
    |                   |  '../Datasets/...0001/MR000005', ...)   |
    | series_reader     | <SimpleITK.SimpleITK.ImageSeriesReader; |
    |                   |  proxy of <Swig Object of type          |
    |                   |  'itk::simple::ImageSeriesReader *' at  |
    |                   |  0x293176760> >                         |
    | image_data        | <SimpleITK.SimpleITK.Image; proxy of    |
    |                   |  <Swig Object of type 'std::vector<     |
    |                   |  itk::simple::Image >::value_type *' at |
    |                   |  0x293089c80> >                         |
    +-------------------+-----------------------------------------+
    +-------------------------+-----------------+
    | Expression              | Result          |
    +-------------------------+-----------------+
    | type(series_ids)        | ⟨class 'tuple'⟩ |
    | len(series_ids)         | 1               |
    | type(series_file_names) | ⟨class 'tuple'⟩ |
    | len(series_file_names)  | 27              |
    | len(image_data)         | 1769472         |
    | image_data.GetSize()    | (256, 256, 27)  |
    +-------------------------+-----------------+
In [23]:
# As shown above, the shape of the original SimpleITK image is (256, 256, 27) instead of the
# common (27, 256, 256) shape, this is due to the different order of the image dimensions,
# so it's necessary to perform the last step to convert the SimpleITK image to a NumPy array

# `GetArrayFromImage` returns a copy of the image data, which can be freely modified
# as it has no effect on the original SimpleITK image
head_mri = sitk.GetArrayFromImage(image_data)

tabulation = Form_Generator()
tabulation.heading_printer("Getting the NumPy array from the SimpleITK image")

statements = ["head_mri = sitk.GetArrayFromImage(image_data)"]
tabulation.statement_generator(statements)

variables = ["head_mri"]
values = [str(reprlib_rules.repr(head_mri))]
tabulation.variable_generator(variables, values)

expressions = [
    "type(head_mri)",
    "len(head_mri)",
    "head_mri.shape",
]
results = [
    str(type(head_mri)),
    str(len(head_mri)),
    str(head_mri.shape),
]
tabulation.expression_generator(expressions, results)
Getting the NumPy array from the SimpleITK image

    +-----------------------------------------------+
    | Statement                                     |
    +-----------------------------------------------+
    | head_mri = sitk.GetArrayFromImage(image_data) |
    +-----------------------------------------------+
    +----------+--------------------------------------------------+
    | Variable | Value                                            |
    +----------+--------------------------------------------------+
    | head_mri | array([[[0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         ...,                                     |
    |          |         [0,...     ...,                          |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0],                 |
    |          |         [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16) |
    +----------+--------------------------------------------------+
    +----------------+-------------------------+
    | Expression     | Result                  |
    +----------------+-------------------------+
    | type(head_mri) | ⟨class 'numpy.ndarray'⟩ |
    | len(head_mri)  | 27                      |
    | head_mri.shape | (27, 256, 256)          |
    +----------------+-------------------------+
In [24]:
head_mri_directory_path = Path(
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
all_paths = sorted(list(head_mri_directory_path.glob("SE*")))
path = all_paths[-1]
new_series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path))
new_series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
    str(path), new_series_ids[0]
)
new_series_reader = sitk.ImageSeriesReader()
new_series_reader.SetFileNames(new_series_file_names)
new_image_data = new_series_reader.Execute()
new_head_mri = sitk.GetArrayFromImage(new_image_data)
# It is easy to see that although the ordered image data can be easily obtained by SimpleITK,
# if other related metadata is needed in addition to the image data, it is still necessary
# to obtain it in the conventional manner
new_all_files = list(path.glob("*"))
new_mri_data = sorted(
    [pydicom.read_file(p) for p in new_all_files],
    key=lambda slice: slice.SliceLocation,
)
new_ordered_slices = {
    index: float(slice.SliceLocation) for index, slice in enumerate(new_mri_data)
}

fig = plt.figure(
    figsize=(figure_size[0], figure_size[1]), constrained_layout=True)

gs = gridspec.GridSpec(
    nrows=2, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[2, 3]
)

ax_1 = fig.add_subplot(gs[0, 0])
grid_DICOM_2D_image(
    head_mri, ax_1, ordered_slices, title="DICOM image grid for MRI axial plane"
)
ax_2 = fig.add_subplot(gs[1, 0])
grid_DICOM_2D_image(
    new_head_mri,
    ax_2,
    new_ordered_slices,
    title="DICOM image grid for MRI coronal plane",
    row_size=8,
    fontsize_adjustment=1,
)

fig.suptitle(
    "Visualization Comparison of Different Plane DICOM Image Grids",
    fontsize="x-large",
    x=0.5,
    y=-0.02,
)

plt.show()
In [25]:
# Create a new function based on the one used earlier, just to get the independent
# plane information
def get_plane(dicom_file):
    IOP = dicom_file.ImageOrientationPatient
    row_xyz = IOP[:3]
    column_xyz = IOP[3:]
    plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
    if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
        return "sagittal"
    elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
        return "coronal"
    elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
        return "axial"
    else:
        "Unknown"


path = all_paths[0]
all_files_0 = list(path.glob("*"))
mri_data_0 = sorted(
    [pydicom.read_file(p) for p in all_files_0],
    key=lambda slice: slice.SliceLocation,
)

fig = plt.figure(
    figsize=(figure_size[0], figure_size[1]), constrained_layout=True)

gs = gridspec.GridSpec(
    nrows=2, ncols=6, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[1, 1]
)

ax_1 = plt.subplot(gs[0, :2])
ax_2 = plt.subplot(gs[0, 2:4])
ax_3 = plt.subplot(gs[0, 4:])
ax_4 = plt.subplot(gs[1, 1:3])
ax_5 = plt.subplot(gs[1, 3:5])
axs = [ax_1, ax_2, ax_3, ax_4, ax_5]

for slice_counter, ax in enumerate(axs):
    DICOM_2D_image(
        ax,
        "bone",
        mri_data_0[slice_counter],
        f"MRI {get_plane(mri_data_0[slice_counter])} plane image at "
        f"{mri_data_0[slice_counter].SliceLocation:^.1f} mm slice location",
        fontsize_adjustment=-1,
        patien_age="52",
    )

fig.suptitle(
    "Visualization Comparison of DICOM Images in Different Planes",
    fontsize="x-large",
    x=0.5,
    y=-0.02,
)

plt.show()

NIfTI in Python¶

Extraction of the NIfTI File¶

In [26]:
path_to_dicom = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
output_folder = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
# The `dicom2nifti.convert_dir` module provides the `convert_directory` function, which is
# automatically imported through the package initialization, to sort all DICOM files one-by-one
# by series, and then convert them to NIfTI format

# By default, the `convert_directory` function is set to True for the `compression` parameter,
# which is used to enable or disable Gzip compression, and for the `reorient` parameter,
# which is used to reorient the DICOM images according to the LAS direction
dicom2nifti.convert_directory(path_to_dicom, output_folder)

# NIfTI stands for Neuroimaging Informatics Technology Initiative and is co-sponsored by
# the National Institute of Mental Health and the National Institute of Neurological Disorders
# and Stroke

# NIfTI defines a neuroimaging data file format designed to meet the needs of the fMRI research
# community, and although some software continues to use other file formats, NIfTI has become
# the most widely used storage standard for fMRI and other MRI research data files

# NIfTI images can be compressed using a standard open-source algorithm called Gzip, which
# greatly reduces file size and therefore the amount of storage required for imaging data

# NIfTI defines a single image file ending in the `.nii` extension, while the compressed
# version appends the `.gz` extension, i.e. `.nii.gz`
nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")

tabulation = Form_Generator()
tabulation.heading_printer(
    "Getting and extracting the NIfTI file through conversion")

statements = [
    """
path_to_dicom = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/SE000001/"
)
output_folder = (
    "../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
    "ST000000/"
)
dicom2nifti.convert_directory(path_to_dicom, output_folder)

nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")
"""
]
tabulation.statement_generator(statements)

variables = ["nifti"]
values = [str(nifti)]
tabulation.variable_generator(variables, values, 3)

expressions = [
    "len(list(nifti.header))",
    "nifti.header['qoffset_x']",
    "nifti.header.get_data_shape()",
    "nifti.shape",
    "nifti.affine",
    "nifti.affine.shape",
]
results = [
    str(len(list(nifti.header))),
    str(nifti.header["qoffset_x"]),
    str(nifti.header.get_data_shape()),
    str(nifti.shape),
    str(nifti.affine),
    str(nifti.affine.shape),
]
tabulation.expression_generator(expressions, results, 3)
Getting and extracting the NIfTI file through conversion

    +-------------------------------------------------------------+
    | Statement                                                   |
    +-------------------------------------------------------------+
    | path_to_dicom = (                                           |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a     |
    |     Normal Male Human Aged 52/"                             |
    |     "ST000000/SE000001/"                                    |
    | )                                                           |
    | output_folder = (                                           |
    |     "../Datasets/An MRI Dicom Data Set of the Head of a     |
    |     Normal Male Human Aged 52/"                             |
    |     "ST000000/"                                             |
    | )                                                           |
    | dicom2nifti.convert_directory(path_to_dicom, output_folder) |
    |                                                             |
    | nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")      |
    +-------------------------------------------------------------+
    +----------+--------------------------------------------------+
    | Variable | Value                                            |
    +----------+--------------------------------------------------+
    | nifti    | ⟨class 'nibabel.nifti1.Nifti1Image'⟩             |
    |          | data shape (256, 256, 27)                        |
    |          | affine:                                          |
    |          | [[-9.36898887e-01  3.20514254e-02                |
    |          |    6.37919828e-02  1.15272324e+02]               |
    |          |  [ 3.03588901e-02  9.27917123e-01                |
    |          |    -8.33337128e-01 -9.72956390e+01]              |
    |          |  [ 1.43172191e-02  1.29802674e-01                |
    |          |    5.94150448e+00 -8.23735046e+01]               |
    |          |  [ 0.00000000e+00  0.00000000e+00                |
    |          |    0.00000000e+00  1.00000000e+00]]              |
    |          | metadata:                                        |
    |          | ⟨class 'nibabel.nifti1.Nifti1Header'⟩ object,    |
    |          |    endian='<'                                    |
    |          | sizeof_hdr      : 348                            |
    |          | data_type       : b''                            |
    |          | db_name         : b''                            |
    |          | extents         : 0                              |
    |          | session_error   : 0                              |
    |          | regular         : b''                            |
    |          | dim_info        : 0                              |
    |          | dim             : [  3 256 256  27   1   1   1   |
    |          |    1]                                            |
    |          | intent_p1       : 0.0                            |
    |          | intent_p2       : 0.0                            |
    |          | intent_p3       : 0.0                            |
    |          | intent_code     : none                           |
    |          | datatype        : uint16                         |
    |          | bitpix          : 16                             |
    |          | slice_start     : 0                              |
    |          | pixdim          : [-1.         0.9375     0.9375 |
    |          |        5.9999995  1.         1.                  |
    |          |   1.         1.       ]                          |
    |          | vox_offset      : 0.0                            |
    |          | scl_slope       : nan                            |
    |          | scl_inter       : nan                            |
    |          | slice_end       : 0                              |
    |          | slice_code      : unknown                        |
    |          | xyzt_units      : 2                              |
    |          | cal_max         : 0.0                            |
    |          | cal_min         : 0.0                            |
    |          | slice_duration  : 0.0                            |
    |          | toffset         : 0.0                            |
    |          | glmax           : 0                              |
    |          | glmin           : 0                              |
    |          | descrip         : b''                            |
    |          | aux_file        : b''                            |
    |          | qform_code      : unknown                        |
    |          | sform_code      : aligned                        |
    |          | quatern_b       : -0.016685799                   |
    |          | quatern_c       : -0.9974202                     |
    |          | quatern_d       : -0.069515765                   |
    |          | qoffset_x       : 115.27232                      |
    |          | qoffset_y       : -97.29564                      |
    |          | qoffset_z       : -82.373505                     |
    |          | srow_x          : [-9.3689889e-01  3.2051425e-02 |
    |          |     6.3791983e-02  1.1527232e+02]                |
    |          | srow_y          : [ 3.035889e-02  9.279171e-01   |
    |          |    -8.333371e-01 -9.729564e+01]                  |
    |          | srow_z          : [ 1.4317219e-02  1.2980267e-01 |
    |          |     5.9415045e+00 -8.2373505e+01]                |
    |          | intent_name     : b''                            |
    |          | magic           : b'n+1'                         |
    +----------+--------------------------------------------------+
    +-------------------------------+---------------------+
    | Expression                    | Result              |
    +-------------------------------+---------------------+
    | len(list(nifti.header))       | 43                  |
    | nifti.header['qoffset_x']     | 115.27232           |
    | nifti.header.get_data_shape() | (256, 256, 27)      |
    | nifti.shape                   | (256, 256, 27)      |
    | nifti.affine                  | [[-9.36898887e-01   |
    |                               |    3.20514254e-02   |
    |                               |    6.37919828e-02   |
    |                               |    1.15272324e+02]  |
    |                               |  [ 3.03588901e-02   |
    |                               |    9.27917123e-01   |
    |                               |    -8.33337128e-01  |
    |                               |    -9.72956390e+01] |
    |                               |  [ 1.43172191e-02   |
    |                               |    1.29802674e-01   |
    |                               |    5.94150448e+00   |
    |                               |    -8.23735046e+01] |
    |                               |  [ 0.00000000e+00   |
    |                               |    0.00000000e+00   |
    |                               |    0.00000000e+00   |
    |                               |    1.00000000e+00]] |
    | nifti.affine.shape            | (4, 4)              |
    +-------------------------------+---------------------+
In [27]:
# Similar to DICOM files, image pixel data can be extracted from the NIfTI file using the
# `get_fdata` function
image_array = nifti.get_fdata()

tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from the NIfTI file")

statements = ["image_array = nifti.get_fdata()"]
tabulation.statement_generator(statements)

variables = ["image_array"]
values = [str(reprlib_rules.repr(image_array))]
tabulation.variable_generator(variables, values)

expressions = [
    "image_array.dtype",
    "image_array.shape",
]
results = [
    str(image_array.dtype),
    str(image_array.shape),
]
tabulation.expression_generator(expressions, results)
Getting pixel data from the NIfTI file

    +---------------------------------+
    | Statement                       |
    +---------------------------------+
    | image_array = nifti.get_fdata() |
    +---------------------------------+
    +-------------+------------------------------------------+
    | Variable    | Value                                    |
    +-------------+------------------------------------------+
    | image_array | array([[[0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |       ... ...,                           |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.],   |
    |             |         [0., 0., 0., ..., 0., 0., 0.]]]) |
    +-------------+------------------------------------------+
    +-------------------+----------------+
    | Expression        | Result         |
    +-------------------+----------------+
    | image_array.dtype | float64        |
    | image_array.shape | (256, 256, 27) |
    +-------------------+----------------+

Visualization of Extracted Images from the NIfTI File¶

In [28]:
def image_display_NIfTI(image, ax, title, rotation=0, cmap="bone"):
    # `scipy.ndimage.rotate` uses spline interpolation of the requested order to return
    # an array rotated in the plane defined by the two axes given by the `axes` parameter

    # By default, the `reshape` parameter is set to True, indicating that the output shape
    # is adjusted so that the input array is fully contained in the output
    rotated_image = ndimage.rotate(image, rotation, reshape=True)
    ax.imshow(rotated_image, cmap)
    ax.grid(False)
    ax.set_title(title, loc="center", pad=10)
    x_ticks = list(range(0, image.shape[1], 50))
    y_ticks = list(range(0, image.shape[0], 50))
    ax.set(xticks=x_ticks, xticklabels=x_ticks,
           yticks=y_ticks, yticklabels=y_ticks)
    ax.set_xlim(left=0)
    ax.set_ylim(top=0)
    ax.minorticks_on()
    return ax


slices_list = np.linspace(
    0, image_array.shape[-1], num=4, endpoint=False, dtype=int)

plt.rcParams["figure.figsize"] = (
    figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)

fig, axs = plt.subplots(nrows=2, ncols=2)

for i in range(4):
    if i == 0:
        title = "Unrotated NIfTI single image"
    else:
        title = f"NIfTI single image rotated {i * 90} degrees"
    image_display_NIfTI(
        image_array[:, :, slices_list[i]],
        axs[i // 2, i % 2],
        title,
        rotation=i * 90,
        cmap="bone",
    )


fig.suptitle(
    "Visual Comparison of NIfTI Single Image at Different Rotation Angles",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [29]:
def NIfTI_title_adder(func):
    @wraps(func)
    def wrapper(*args, rotation=0, title=None, **kwargs):
        if title is None:
            if rotation == 0:
                title = "Unrotated NIfTI image grid"
            else:
                title = f"NIfTI image grid with rotated {rotation} degrees"
        ax = func(*args, rotation, **kwargs)
        ax.set_title(
            title,
            loc="center",
            pad=10,
        )

    return wrapper


@NIfTI_title_adder
def grid_NIfTI_2D_image(
    inputs,
    ax,
    rotation=0,
    slice_axis=0,
    slice_interval=1,
    n=None,
    row_size=9,
    pad_value=1,
):
    # As in the case of the SimpleITK package, NIfTI data converted from DICOM data changes
    # the dimensional order of the 3D volume after conversion, resulting in a shape of
    # (256, 256, 27) instead of the common (27, 256, 256) shape

    # In fact, the 3D image in NIfTI data (3D or 4D) can be regarded as a stack of 2D images
    # along three axes respectively, so it seems that the order of the three dimensions, i.e.,
    # the three axes, is no longer so important here, but the important thing is the way to
    # extract the stacked 2D images along any axis
    if slice_axis == 2:
        inputs = inputs.transpose(2, 0, 1)
    elif slice_axis == 1:
        inputs = inputs.transpose(1, 0, 2)
    if n is None:
        n = math.ceil(inputs.shape[0] / slice_interval)
    else:
        n = math.ceil(n / slice_interval)
    nrows = math.ceil(n / row_size)
    cmap = "bone"

    rotated_inputs = [
        ndimage.rotate(inputs[i, :, :], rotation, reshape=True)
        for i in range(len(inputs))
    ]
    rotated_inputs = np.array(rotated_inputs)
    inputs = torch.Tensor(rotated_inputs.astype(int)).unsqueeze(1)
    inputs = inputs[::slice_interval, :, :, :]

    for i in range(nrows):
        if len(inputs) > row_size:
            row_images = inputs[:row_size]
            inputs = inputs[row_size:]
        else:
            row_images = inputs[:]
        grid_row = make_grid(
            row_images, nrow=row_size, normalize=True, pad_value=pad_value
        )
        scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
        if i == 0:
            if nrows != 1:
                scalar_image = scalar_row
            else:
                scalar_image = np.full(
                    (
                        scalar_row.shape[0],
                        math.ceil(
                            scalar_row.shape[1] / inputs.shape[0] * row_size),
                    ),
                    np.nan,
                )
                scalar_image[:, : scalar_row.shape[1]] = scalar_row
        elif scalar_image.shape[1] == scalar_row.shape[1]:
            scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
        else:
            image_to_concat = np.full(
                (scalar_image.shape[0] +
                 scalar_row.shape[0], scalar_image.shape[1]),
                np.nan,
            )
            image_to_concat[: scalar_image.shape[0], :] = scalar_image
            image_to_concat[scalar_image.shape[0]:,
                            : scalar_row.shape[1]] = scalar_row
            scalar_image = image_to_concat
    ax.imshow(scalar_image, cmap=cmap)
    ax.set(xticks=[], yticks=[], frame_on=False)
    ax.grid(False)
    return ax


fig, axs = plt.subplots(4, 1, figsize=(
    figure_size[0], figure_size[1] / 7 * 12))

for i in range(4):
    grid_NIfTI_2D_image(image_array, axs[i], rotation=i * 90, slice_axis=2)


fig.suptitle(
    "Visual Comparison of NIfTI Image Grid at Different Rotation Angles",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()

Writing of the NIfTI file¶

In [30]:
# In this preprocessing step, a simple threshold is set to turn all image voxels with values
# not greater than 300 to 0

# The principle is to determine the True and False values of all image voxels and multiply them
# by the corresponding value of 1 or 0, respectively
image_array_processed = image_array * (image_array > 300)

tabulation = Form_Generator()
tabulation.heading_printer("A simple preprocessing step")

statements = ["image_array_processed = image_array * (image_array > 300)"]
tabulation.statement_generator(statements)

variables = ["image_array_processed"]
values = [str(reprlib_rules.repr(image_array_processed))]
tabulation.variable_generator(variables, values, 9)

expressions = [
    "image_array > 300",
    "image_array.dtype",
    "image_array.shape",
]
results = [
    str(reprlib_rules.repr(image_array > 300)),
    str(image_array_processed.dtype),
    str(image_array_processed.shape),
]
tabulation.expression_generator(expressions, results, 9)
A simple preprocessing step

    +-----------------------------------------------------------+
    | Statement                                                 |
    +-----------------------------------------------------------+
    | image_array_processed = image_array * (image_array > 300) |
    +-----------------------------------------------------------+
    +-----------------------+-----------------------------------+
    | Variable              | Value                             |
    +-----------------------+-----------------------------------+
    | image_array_processed | array([[[0., 0., 0., ..., 0., 0., |
    |                       |          0.],                     |
    |                       |         [0., 0., 0., ..., 0., 0., |
    |                       |          0.],                     |
    |                       |         [0., 0., 0., ..., 0., 0., |
    |                       |          0.],                     |
    |                       |       ... ...,                    |
    |                       |         [0., 0., 0., ..., 0., 0., |
    |                       |          0.],                     |
    |                       |         [0., 0., 0., ..., 0., 0., |
    |                       |          0.],                     |
    |                       |         [0., 0., 0., ..., 0., 0., |
    |                       |          0.]]])                   |
    +-----------------------+-----------------------------------+
    +-------------------+------------------------------------+
    | Expression        | Result                             |
    +-------------------+------------------------------------+
    | image_array > 300 | array([[[False, False, False, ..., |
    |                   |          False, False, False],     |
    |                   |         [False, False, False, ..., |
    |                   |          False, False, False],     |
    |                   |         [... False],               |
    |                   |         [False, False, False, ..., |
    |                   |          False, False, False],     |
    |                   |         [False, False, False, ..., |
    |                   |          False, False, False]]])   |
    | image_array.dtype | float64                            |
    | image_array.shape | (256, 256, 27)                     |
    +-------------------+------------------------------------+
In [31]:
plt.rcParams["figure.figsize"] = (
    figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)

fig, axs = plt.subplots(nrows=2, ncols=2)

for i in range(4):
    if i == 0:
        title = "Unrotated preprocessed NIfTI single image"
    else:
        title = f"Preprocessed NIfTI single image rotated {i * 90} degrees"
    image_display_NIfTI(
        image_array_processed[:, :, slices_list[i]],
        axs[i // 2, i % 2],
        title,
        rotation=i * 90,
        cmap="bone",
    )


fig.suptitle(
    "Visual Comparison of Preprocessed NIfTI Single Images with Different Rotation Angles",
    fontsize="x-large",
    x=0.5,
    y=0,
)

plt.tight_layout()
plt.show()
In [32]:
# `nib.Nifti1Image` returns the initialized image, which is a combination of array-like object,
# affine matrix, header and optional metadata in `extra`, et cetera

# A NiBabel image consists of three parts: a 3D or 4D array of image data, an affine matrix
# that tells the position of the image array data in a reference space, and image metadata
# to describe the image

# In particular, the affine matrix is an identity affine that stores the relationship between
# voxel coordinates in the image data array and coordinates in reference space, where the
# voxel axes are aligned with the world axes
processed_nifti = nib.Nifti1Image(image_array_processed, nifti.affine)

# `nib.save` saves NIfTI image with user-defined filename
nib.save(processed_nifti, f"{output_folder}201_t2w_tse_processed.nii.gz")

tabulation = Form_Generator()
tabulation.heading_printer("Initialization and saving of NIfTI image")

statements = [
    """
processed_nifti = nib.Nifti1Image(image_array_processed, nifti.affine)

nib.save(processed_nifti, f"{output_folder}201_t2w_tse_processed.nii.gz")
"""
]
tabulation.statement_generator(statements)

variables = ["processed_nifti"]
values = [str(processed_nifti)]
tabulation.variable_generator(variables, values, 3)

expressions = [
    "len(list(processed_nifti.header))",
    "processed_nifti.shape",
    "processed_nifti.affine",
    "processed_nifti.affine.shape",
    "processed_nifti.extra",
]
results = [
    str(len(list(processed_nifti.header))),
    str(processed_nifti.shape),
    str(processed_nifti.affine),
    str(processed_nifti.affine.shape),
    str(processed_nifti.extra),
]
tabulation.expression_generator(expressions, results, 3)
Initialization and saving of NIfTI image

    +----------------------------------------------------------+
    | Statement                                                |
    +----------------------------------------------------------+
    | processed_nifti = nib.Nifti1Image(image_array_processed, |
    |     nifti.affine)                                        |
    |                                                          |
    | nib.save(processed_nifti,                                |
    |     f"{output_folder}201_t2w_tse_processed.nii.gz")      |
    +----------------------------------------------------------+
    +-----------------+-------------------------------------------+
    | Variable        | Value                                     |
    +-----------------+-------------------------------------------+
    | processed_nifti | ⟨class 'nibabel.nifti1.Nifti1Image'⟩      |
    |                 | data shape (256, 256, 27)                 |
    |                 | affine:                                   |
    |                 | [[-9.36898887e-01  3.20514254e-02         |
    |                 |    6.37919828e-02  1.15272324e+02]        |
    |                 |  [ 3.03588901e-02  9.27917123e-01         |
    |                 |    -8.33337128e-01 -9.72956390e+01]       |
    |                 |  [ 1.43172191e-02  1.29802674e-01         |
    |                 |    5.94150448e+00 -8.23735046e+01]        |
    |                 |  [ 0.00000000e+00  0.00000000e+00         |
    |                 |    0.00000000e+00  1.00000000e+00]]       |
    |                 | metadata:                                 |
    |                 | ⟨class 'nibabel.nifti1.Nifti1Header'⟩     |
    |                 |    object, endian='<'                     |
    |                 | sizeof_hdr      : 348                     |
    |                 | data_type       : b''                     |
    |                 | db_name         : b''                     |
    |                 | extents         : 0                       |
    |                 | session_error   : 0                       |
    |                 | regular         : b''                     |
    |                 | dim_info        : 0                       |
    |                 | dim             : [  3 256 256  27   1    |
    |                 |    1   1   1]                             |
    |                 | intent_p1       : 0.0                     |
    |                 | intent_p2       : 0.0                     |
    |                 | intent_p3       : 0.0                     |
    |                 | intent_code     : none                    |
    |                 | datatype        : float64                 |
    |                 | bitpix          : 64                      |
    |                 | slice_start     : 0                       |
    |                 | pixdim          : [-1.                    |
    |                 |    0.93749994  0.9375      5.9999995   1. |
    |                 |             1.                            |
    |                 |   1.          1.        ]                 |
    |                 | vox_offset      : 0.0                     |
    |                 | scl_slope       : nan                     |
    |                 | scl_inter       : nan                     |
    |                 | slice_end       : 0                       |
    |                 | slice_code      : unknown                 |
    |                 | xyzt_units      : 0                       |
    |                 | cal_max         : 0.0                     |
    |                 | cal_min         : 0.0                     |
    |                 | slice_duration  : 0.0                     |
    |                 | toffset         : 0.0                     |
    |                 | glmax           : 0                       |
    |                 | glmin           : 0                       |
    |                 | descrip         : b''                     |
    |                 | aux_file        : b''                     |
    |                 | qform_code      : unknown                 |
    |                 | sform_code      : aligned                 |
    |                 | quatern_b       : -0.016685799            |
    |                 | quatern_c       : -0.9974202              |
    |                 | quatern_d       : -0.06951577             |
    |                 | qoffset_x       : 115.27232               |
    |                 | qoffset_y       : -97.29564               |
    |                 | qoffset_z       : -82.373505              |
    |                 | srow_x          : [-9.3689889e-01         |
    |                 |    3.2051425e-02  6.3791983e-02           |
    |                 |    1.1527232e+02]                         |
    |                 | srow_y          : [ 3.035889e-02          |
    |                 |    9.279171e-01 -8.333371e-01             |
    |                 |    -9.729564e+01]                         |
    |                 | srow_z          : [ 1.4317219e-02         |
    |                 |    1.2980267e-01  5.9415045e+00           |
    |                 |    -8.2373505e+01]                        |
    |                 | intent_name     : b''                     |
    |                 | magic           : b'n+1'                  |
    +-----------------+-------------------------------------------+
    +-----------------------------------+---------------------+
    | Expression                        | Result              |
    +-----------------------------------+---------------------+
    | len(list(processed_nifti.header)) | 43                  |
    | processed_nifti.shape             | (256, 256, 27)      |
    | processed_nifti.affine            | [[-9.36898887e-01   |
    |                                   |    3.20514254e-02   |
    |                                   |    6.37919828e-02   |
    |                                   |    1.15272324e+02]  |
    |                                   |  [ 3.03588901e-02   |
    |                                   |    9.27917123e-01   |
    |                                   |    -8.33337128e-01  |
    |                                   |    -9.72956390e+01] |
    |                                   |  [ 1.43172191e-02   |
    |                                   |    1.29802674e-01   |
    |                                   |    5.94150448e+00   |
    |                                   |    -8.23735046e+01] |
    |                                   |  [ 0.00000000e+00   |
    |                                   |    0.00000000e+00   |
    |                                   |    0.00000000e+00   |
    |                                   |    1.00000000e+00]] |
    | processed_nifti.affine.shape      | (4, 4)              |
    | processed_nifti.extra             | {}                  |
    +-----------------------------------+---------------------+

Preprocessing in Python¶

Orientation¶

In [33]:
# The dataset used in this practice project was collected as part of the Information eXtraction
# from Images (IXI) project, which collects MRI data from nearly 600 normal healthy subjects

# The MRI acquisition protocol for each subject consisted of T1, T2, and PD-weighted images,
# MRA images, and diffusion-weighted images (15 directions)

# The data in this dataset were acquired at three different hospitals in London using 1.5T and
# 3T MRI scanning systems

# This dataset was provided by IXI Dataset (https://brain-development.org/ixi-dataset/),
# license type is Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0)
# https://creativecommons.org/licenses/by-sa/3.0/legalcode
brain_mri_path = "../Datasets/IXI-T1/"
brain_mri = nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz")
brain_mri_data = brain_mri.get_fdata()
shape = brain_mri.shape
# The affine matrix A describes the mapping from pixel coordinates to scanner or world
# coordinates and is a 4 x 4 matrix that can be represented as follows:
# A = | m_11, m_12, m_13, a |
#     | m_21, m_22, m_23, b |
#     | m_31, m_32, m_33, c |
#     |    0,    0,    0, 1 |

# In the affine matrix A, a 3 x 3 matrix of m-series elements is responsible for rotating,
# scaling, and shearing the coordinates, denoted by M

# For a detailed description of rotation, scaling, translation, or shearing, see the following
# links: https://en.wikipedia.org/wiki/Affine_transformation#Image_transformation

# In the affine matrix A, there is also a 3 x 1 translation vector (a, b, c) that is
# responsible for the translation or offset of the affine matrix in three dimensions

# In this way, the pixel coordinates (i, j, k) are mapped to the coordinates (x, y, z) of the
# scanner or the world as follows:
# | x | = | m_11, m_12, m_13 | * | i | + | a' | = M * | i | + | a' |
# | y |   | m_21, m_22, m_23 |   | j |   | b' |       | j |   | b' |
# | z |   | m_31, m_32, m_33 |   | k |   | c' |       | k |   | c' |,
# where a' = a * i, b' = b * j, c' = c * k

# The addition of the extra row [0, 0, 0, 1] is a trick to put the translation part into the
# same matrix A as the matrix M, so that both parts can be applied by matrix multiplication

# To implement the above trick, an extra 1 needs to be added to the input (pixel) and output
# (scanner or world) coordinate vectors, like this:
# | x | = | m_11, m_12, m_13, a | * | i | = A * | i |
# | y |   | m_21, m_22, m_23, b |   | j |       | j |
# | z |   | m_31, m_32, m_33, c |   | k |       | k |
# | 1 |   |    0,    0,    0, 1 |   | 1 |       | 1 |

# A more detailed explanation can be found at this link:
# https://nipy.org/nibabel/coordinate_systems.html#voxel-coordinates-are-in-voxel-space
affine = brain_mri.affine

tabulation = Form_Generator()
tabulation.heading_printer(
    "Loading MRI data in NIfTI format and obtaining its affine matrix"
)

statements = [
    """
brain_mri_path = "../Datasets/IXI-T1/"
brain_mri = nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz")
brain_mri_data = brain_mri.get_fdata()
shape = brain_mri.shape
affine = brain_mri.affine
"""
]
tabulation.statement_generator(statements)

variables = ["brain_mri", "brain_mri_data", "shape", "affine"]
values = [
    str(brain_mri),
    str(reprlib_rules.repr(brain_mri_data)),
    str(shape),
    str(affine),
]
tabulation.variable_generator(variables, values, 3)

expressions = [
    "len(list(brain_mri.header))",
    "brain_mri.header.get_zooms()",
    "affine.shape",
    "nib.aff2axcodes(affine)",
    "affine[:3, :3]",
    "affine[:3, 3:]",
    "affine[3:, :]",
]

results = [
    str(len(list(brain_mri.header))),
    str(brain_mri.header.get_zooms()),
    str(affine.shape),
    str(nib.aff2axcodes(affine)),
    str(affine[:3, :3]),
    str(affine[:3, 3:]),
    str(affine[3:, :]),
]
tabulation.expression_generator(expressions, results, 3)
Loading MRI data in NIfTI format and obtaining its affine matrix

    +-------------------------------------------------------------+
    | Statement                                                   |
    +-------------------------------------------------------------+
    | brain_mri_path = "../Datasets/IXI-T1/"                      |
    | brain_mri =                                                 |
    |     nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz") |
    | brain_mri_data = brain_mri.get_fdata()                      |
    | shape = brain_mri.shape                                     |
    | affine = brain_mri.affine                                   |
    +-------------------------------------------------------------+
    +----------------+--------------------------------------------+
    | Variable       | Value                                      |
    +----------------+--------------------------------------------+
    | brain_mri      | ⟨class 'nibabel.nifti1.Nifti1Image'⟩       |
    |                | data shape (256, 256, 150)                 |
    |                | affine:                                    |
    |                | [[ 1.89821944e-02 -2.72075552e-03          |
    |                |    1.19975281e+00 -9.06798553e+01]         |
    |                |  [-9.27821696e-01  1.32986516e-01          |
    |                |    2.45456006e-02  1.02829445e+02]         |
    |                |  [ 1.33014351e-01  9.28015888e-01          |
    |                |    5.71511449e-11 -1.14823784e+02]         |
    |                |  [ 0.00000000e+00  0.00000000e+00          |
    |                |    0.00000000e+00  1.00000000e+00]]        |
    |                | metadata:                                  |
    |                | ⟨class 'nibabel.nifti1.Nifti1Header'⟩      |
    |                |    object, endian='<'                      |
    |                | sizeof_hdr      : 348                      |
    |                | data_type       : b''                      |
    |                | db_name         : b''                      |
    |                | extents         : 0                        |
    |                | session_error   : 0                        |
    |                | regular         : b'r'                     |
    |                | dim_info        : 0                        |
    |                | dim             : [  3 256 256 150   1   1 |
    |                |      1   1]                                |
    |                | intent_p1       : 0.0                      |
    |                | intent_p2       : 0.0                      |
    |                | intent_p3       : 0.0                      |
    |                | intent_code     : none                     |
    |                | datatype        : int16                    |
    |                | bitpix          : 16                       |
    |                | slice_start     : 0                        |
    |                | pixdim          : [-1.         0.9375      |
    |                |    0.9375     1.2000039  0.         0.     |
    |                |   0.         0.       ]                    |
    |                | vox_offset      : 0.0                      |
    |                | scl_slope       : nan                      |
    |                | scl_inter       : nan                      |
    |                | slice_end       : 0                        |
    |                | slice_code      : unknown                  |
    |                | xyzt_units      : 10                       |
    |                | cal_max         : 0.0                      |
    |                | cal_min         : 0.0                      |
    |                | slice_duration  : 0.0                      |
    |                | toffset         : 0.0                      |
    |                | glmax           : 0                        |
    |                | glmin           : 0                        |
    |                | descrip         : b'MR'                    |
    |                | aux_file        : b''                      |
    |                | qform_code      : scanner                  |
    |                | sform_code      : scanner                  |
    |                | quatern_b       : 0.46861374               |
    |                | quatern_c       : -0.52952915              |
    |                | quatern_d       : -0.4576844               |
    |                | qoffset_x       : -90.679855               |
    |                | qoffset_y       : 102.829445               |
    |                | qoffset_z       : -114.823784              |
    |                | srow_x          : [ 1.8982194e-02          |
    |                |    -2.7207555e-03  1.1997528e+00           |
    |                |    -9.0679855e+01]                         |
    |                | srow_y          : [-9.27821696e-01         |
    |                |    1.32986516e-01  2.45456006e-02          |
    |                |    1.02829445e+02]                         |
    |                | srow_z          : [ 1.33014351e-01         |
    |                |    9.28015888e-01  5.71511449e-11          |
    |                |    -1.14823784e+02]                        |
    |                | intent_name     : b''                      |
    |                | magic           : b'n+1'                   |
    | brain_mri_data | array([[[0., 0., 0., ..., 0., 0., 0.],     |
    |                |         [0., 0., 0., ..., 0., 0., 0.],     |
    |                |         [0., 0., 0., ..., 0., 0., 0.],     |
    |                |       ... ...,                             |
    |                |         [0., 0., 0., ..., 0., 0., 0.],     |
    |                |         [0., 0., 0., ..., 0., 0., 0.],     |
    |                |         [0., 0., 0., ..., 0., 0., 0.]]])   |
    | shape          | (256, 256, 150)                            |
    | affine         | [[ 1.89821944e-02 -2.72075552e-03          |
    |                |    1.19975281e+00 -9.06798553e+01]         |
    |                |  [-9.27821696e-01  1.32986516e-01          |
    |                |    2.45456006e-02  1.02829445e+02]         |
    |                |  [ 1.33014351e-01  9.28015888e-01          |
    |                |    5.71511449e-11 -1.14823784e+02]         |
    |                |  [ 0.00000000e+00  0.00000000e+00          |
    |                |    0.00000000e+00  1.00000000e+00]]        |
    +----------------+--------------------------------------------+
    +------------------------------+-----------------------------+
    | Expression                   | Result                      |
    +------------------------------+-----------------------------+
    | len(list(brain_mri.header))  | 43                          |
    | brain_mri.header.get_zooms() | (0.9375, 0.9375, 1.2000039) |
    | affine.shape                 | (4, 4)                      |
    | nib.aff2axcodes(affine)      | ('P', 'S', 'R')             |
    | affine[:3, :3]               | [[ 1.89821944e-02           |
    |                              |    -2.72075552e-03          |
    |                              |    1.19975281e+00]          |
    |                              |  [-9.27821696e-01           |
    |                              |    1.32986516e-01           |
    |                              |    2.45456006e-02]          |
    |                              |  [ 1.33014351e-01           |
    |                              |    9.28015888e-01           |
    |                              |    5.71511449e-11]]         |
    | affine[:3, 3:]               | [[ -90.67985535]            |
    |                              |  [ 102.82944489]            |
    |                              |  [-114.82378387]]           |
    | affine[3:, :]                | [[0. 0. 0. 1.]]             |
    +------------------------------+-----------------------------+
In [34]:
fig = plt.figure(
    figsize=(figure_size[0], figure_size[1] / 7 * 16), constrained_layout=True
)

gs = gridspec.GridSpec(
    nrows=3, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[50, 50, 49]
)

ax_1 = plt.subplot(gs[0, 0])
ax_2 = plt.subplot(gs[1, 0])
ax_3 = plt.subplot(gs[2, 0])
axs = [ax_1, ax_2, ax_3]

for i in range(3):
    # The image voxel orientation indicates how the image array data is stored and is
    # represented by three axes: left/right, anterior/posterior, and inferior/superior

    # The image orientation can be obtained from the affine matrix of NIfTI data, which is
    # a tuple of three single-letter strings representing the order and increasing direction
    # of the axes

    # For example, the orientation indicated by ('P', 'S', 'R') corresponds to the first axis
    # ordered from anterior to posterior (coronal plane), the second axis ordered from inferior
    # to superior (axial plane), and the third axis ordered from left to right (sagittal plane)
    plane_dict = {
        "S": "axial",
        "I": "axial",
        "L": "sagittal",
        "R": "sagittal",
        "A": "coronal",
        "P": "coronal",
    }
    nth = {1: "first", 2: "second", 3: "third"}
    # `nib.aff2axcodes` returns the axis orientation codes for the affine matrix, expressed
    # as a tuple of labels for the positive end of voxel axes
    title = f"NIfTI image grid on the {plane_dict[nib.aff2axcodes(affine)[i]]} axis (the {nth[i+1]} axis)"
    slice_interval, row_size = (8, 8) if shape[i] == 256 else (5, 6)
    # The 3D image can be viewed as a stack of planes on three different axis planes, each axis
    # plane corresponding to a different pairwise combination of the three coordinate axes

    # Since the first two axes have a common coordinate axis in their respective paired
    # coordinate axes, their planar orientations are imaged identically, and the planar image
    # of the third axis will be rotated relative to the planes of the first two axes

    # Although the custom function provides a way to visualize rotated sub-images, it is also
    # possible to obtain sub-images in the target direction by simply swapping the axes and
    # flipping the axis ordering, but be aware that neither of these methods change the
    # affine matrix of the original data, so it should be used with caution and only for
    # visualization purposes
    if i == 0:
        # `np.flip` reverses the order of the elements in the array along the given axes, where
        # the argument `axis` represents the axis to be flipped, and the default value is None,
        # which means that all axes of the input array will be flipped
        new_brain_mri = np.flip(brain_mri_data, axis=1)
        new_brain_mri = np.flip(new_brain_mri, axis=2)
    elif i == 1:
        new_brain_mri = np.flip(brain_mri_data, 0)
        new_brain_mri = np.flip(new_brain_mri, 2)
    else:
        # `np.swapaxes` interchanges the two axes of the array, where the arguments `axis1` and
        # `axis2` represent the first axis and the second axis, respectively
        brain_mri_swapped = np.swapaxes(brain_mri_data, axis1=0, axis2=1)
        new_brain_mri = np.flip(brain_mri_swapped, 0)
    grid_NIfTI_2D_image(
        new_brain_mri,
        axs[i],
        slice_axis=i,
        slice_interval=slice_interval,
        row_size=row_size,
        title=title,
    )


fig.suptitle(
    "Visual Comparison of NIfTI Image Grids on Three Different Axes",
    fontsize="x-large",
    x=0.5,
    y=-0.01,
)

plt.show()
In [35]:
# To calculate the physical coordinates of the offset, i.e., the position of the voxel coordinates (0, 0, 0)
# in physical space, it is sufficient to multiply the affine matrix by these coordinates, but note that
# for convenience it is necessary to append a 1 when using a 4 x 4 matrix
voxel_coord = np.array((0, 0, 0, 1))
# The `@` operator is used for matrix multiplication and is equivalent to the `np. matmul`
# function that handles the multiplication of two arrays
physical_coord = affine @ voxel_coord

# It is also possible to skip the shortcut and calculate the result by hand, in which case
# there is no need to append 1, since the core matrix in the calculation is a 3 x 3 matrix
voxel_coord_manual = np.array((0, 0, 0))
physical_coord_manual = affine[:3, :3] @ voxel_coord_manual
physical_coord_manual = physical_coord_manual + affine[:3, 3]

# If it is necessary to convert physical coordinates to the pixel/voxel coordinates,
# the inverse of the affine matrix needs to be computed using `np.linalg.inv`, and then
# this inverse is multiplied by the physical coordinates

# `np.linalg.inv` computes the (multiplicative) inverse of a matrix
voxel_coords = np.linalg.inv(affine) @ physical_coord

tabulation = Form_Generator()
tabulation.heading_printer(
    "Calculation of coordinate transformations between coordinate systems"
)

statements = [
    """
voxel_coord = np.array((0, 0, 0, 1))
physical_coord = affine @ voxel_coord

voxel_coord_manual = np.array((0, 0, 0))
physical_coord_manual = affine[:3, :3] @ voxel_coord_manual
physical_coord_manual = physical_coord_manual + affine[:3, 3]

voxel_coords = np.linalg.inv(affine) @ physical_coord
"""
]
tabulation.statement_generator(statements)

variables = [
    "voxel_coord",
    "physical_coord",
    "voxel_coord_manual",
    "physical_coord_manual",
    "voxel_coords",
]
values = [
    str(voxel_coord),
    str(physical_coord),
    str(voxel_coord_manual),
    str(physical_coord_manual),
    str(voxel_coords),
]
tabulation.variable_generator(variables, values, 1)

expressions = [
    "voxel_coord.shape",
    "physical_coord.shape",
    "affine[:, 3:]",
    "affine[:, 3:].shape",
    "voxel_coord_manual.shape",
    "physical_coord_manual.shape",
    "voxel_coords.round()",
    "voxel_coords.shape",
]

results = [
    str(voxel_coord.shape),
    str(physical_coord.shape),
    str(affine[:, 3:]),
    str(affine[:, 3:].shape),
    str(voxel_coord_manual.shape),
    str(physical_coord_manual.shape),
    str(voxel_coords.round()),
    str(voxel_coords.shape),
]
tabulation.expression_generator(expressions, results, 1)
Calculation of coordinate transformations between coordinate systems

    +-------------------------------------------------------------+
    | Statement                                                   |
    +-------------------------------------------------------------+
    | voxel_coord = np.array((0, 0, 0, 1))                        |
    | physical_coord = affine @ voxel_coord                       |
    |                                                             |
    | voxel_coord_manual = np.array((0, 0, 0))                    |
    | physical_coord_manual = affine[:3, :3] @ voxel_coord_manual |
    | physical_coord_manual = physical_coord_manual + affine[:3,  |
    |     3]                                                      |
    |                                                             |
    | voxel_coords = np.linalg.inv(affine) @ physical_coord       |
    +-------------------------------------------------------------+
    +-----------------------+-----------------------------------+
    | Variable              | Value                             |
    +-----------------------+-----------------------------------+
    | voxel_coord           | [0 0 0 1]                         |
    | physical_coord        | [ -90.67985535  102.82944489      |
    |                       |  -114.82378387    1.        ]     |
    | voxel_coord_manual    | [0 0 0]                           |
    | physical_coord_manual | [ -90.67985535  102.82944489      |
    |                       |  -114.82378387]                   |
    | voxel_coords          | [-7.10542736e-15  0.00000000e+00  |
    |                       |  0.00000000e+00  1.00000000e+00]  |
    +-----------------------+-----------------------------------+
    +-----------------------------+-------------------+
    | Expression                  | Result            |
    +-----------------------------+-------------------+
    | voxel_coord.shape           | (4,)              |
    | physical_coord.shape        | (4,)              |
    | affine[:, 3:]               | [[ -90.67985535]  |
    |                             |  [ 102.82944489]  |
    |                             |  [-114.82378387]  |
    |                             |  [   1.        ]] |
    | affine[:, 3:].shape         | (4, 1)            |
    | voxel_coord_manual.shape    | (3,)              |
    | physical_coord_manual.shape | (3,)              |
    | voxel_coords.round()        | [-0.  0.  0.  1.] |
    | voxel_coords.shape          | (4,)              |
    +-----------------------------+-------------------+

Canonicalization¶

In [36]:
# By convention, the world axes of NiBabel are always in the RAS+ orientation, i.e.,
# left to right, posterior to anterior, and inferior to superior, and this voxel orientation
# is referred to as the canonical voxel orientation

# Sometimes it may be necessary to rearrange the image voxel axes to be as close as possible to
# the RAS+ orientation, since RAS+ is the canonical world orientation

# Note that rearranging the voxel axes means reversing and/or rearranging the voxel axes

# `nib.as_closest_canonical` returns the closest canonical image after reordering the data
brain_mri_canonical = nib.as_closest_canonical(brain_mri)
brain_mri_canonical_data = brain_mri_canonical.get_fdata()
canonical_affine = brain_mri_canonical.affine

tabulation = Form_Generator()
tabulation.heading_printer("Canonicalization of NIfTI data")

statements = [
    """
brain_mri_canonical = nib.as_closest_canonical(brain_mri)
brain_mri_canonical_data = brain_mri_canonical.get_fdata()
canonical_affine = brain_mri_canonical.affine
"""
]
tabulation.statement_generator(statements)

variables = ["brain_mri_canonical", "brain_mri_canonical_data", "canonical_affine"]
values = [
    str(brain_mri_canonical),
    str(reprlib_rules.repr(brain_mri_canonical_data)),
    str(canonical_affine),
]
tabulation.variable_generator(variables, values, 2)

expressions = ["brain_mri_canonical_data.shape", "nib.aff2axcodes(canonical_affine)"]

results = [str(brain_mri_canonical_data.shape), str(nib.aff2axcodes(canonical_affine))]
tabulation.expression_generator(expressions, results)
Canonicalization of NIfTI data

    +------------------------------------------------------------+
    | Statement                                                  |
    +------------------------------------------------------------+
    | brain_mri_canonical = nib.as_closest_canonical(brain_mri)  |
    | brain_mri_canonical_data = brain_mri_canonical.get_fdata() |
    | canonical_affine = brain_mri_canonical.affine              |
    +------------------------------------------------------------+
    +--------------------------+----------------------------------+
    | Variable                 | Value                            |
    +--------------------------+----------------------------------+
    | brain_mri_canonical      | <class                           |
    |                          |   'nibabel.nifti1.Nifti1Image'>  |
    |                          | data shape (150, 256, 256)       |
    |                          | affine:                          |
    |                          | [[ 1.19975281e+00                |
    |                          |   -1.89821944e-02                |
    |                          |   -2.72075552e-03                |
    |                          |   -8.58393958e+01]               |
    |                          |  [ 2.45456006e-02                |
    |                          |   9.27821696e-01  1.32986516e-01 |
    |                          |   -1.33765088e+02]               |
    |                          |  [ 5.71511449e-11                |
    |                          |   -1.33014351e-01                |
    |                          |   9.28015888e-01                 |
    |                          |   -8.09051243e+01]               |
    |                          |  [ 0.00000000e+00                |
    |                          |   0.00000000e+00  0.00000000e+00 |
    |                          |    1.00000000e+00]]              |
    |                          | metadata:                        |
    |                          | <class                           |
    |                          |   'nibabel.nifti1.Nifti1Header'> |
    |                          |   object, endian='<'             |
    |                          | sizeof_hdr      : 348            |
    |                          | data_type       : b''            |
    |                          | db_name         : b''            |
    |                          | extents         : 0              |
    |                          | session_error   : 0              |
    |                          | regular         : b'r'           |
    |                          | dim_info        : 0              |
    |                          | dim             : [  3 150 256   |
    |                          |   256   1   1   1   1]           |
    |                          | intent_p1       : 0.0            |
    |                          | intent_p2       : 0.0            |
    |                          | intent_p3       : 0.0            |
    |                          | intent_code     : none           |
    |                          | datatype        : int16          |
    |                          | bitpix          : 16             |
    |                          | slice_start     : 0              |
    |                          | pixdim          : [1.            |
    |                          |   1.2000039 0.9375    0.9375     |
    |                          |   1.        1.        1.         |
    |                          |  1.       ]                      |
    |                          | vox_offset      : 0.0            |
    |                          | scl_slope       : nan            |
    |                          | scl_inter       : nan            |
    |                          | slice_end       : 0              |
    |                          | slice_code      : unknown        |
    |                          | xyzt_units      : 10             |
    |                          | cal_max         : 0.0            |
    |                          | cal_min         : 0.0            |
    |                          | slice_duration  : 0.0            |
    |                          | toffset         : 0.0            |
    |                          | glmax           : 0              |
    |                          | glmin           : 0              |
    |                          | descrip         : b'MR'          |
    |                          | aux_file        : b''            |
    |                          | qform_code      : unknown        |
    |                          | sform_code      : aligned        |
    |                          | quatern_b       : -0.071117364   |
    |                          | quatern_c       : -0.0007274148  |
    |                          | quatern_d       : 0.010201936    |
    |                          | qoffset_x       : -85.83939      |
    |                          | qoffset_y       : -133.76509     |
    |                          | qoffset_z       : -80.90512      |
    |                          | srow_x          : [              |
    |                          |   1.1997528e+00 -1.8982194e-02   |
    |                          |   -2.7207555e-03 -8.5839394e+01] |
    |                          | srow_y          : [              |
    |                          |   2.4545601e-02  9.2782170e-01   |
    |                          |   1.3298652e-01 -1.3376509e+02]  |
    |                          | srow_z          : [              |
    |                          |   5.7151145e-11 -1.3301435e-01   |
    |                          |   9.2801589e-01 -8.0905121e+01]  |
    |                          | intent_name     : b''            |
    |                          | magic           : b'n+1'         |
    | brain_mri_canonical_data | array([[[0., 0., 0., ..., 0.,    |
    |                          |   0., 0.],                       |
    |                          |         [0., 0., 0., ..., 0.,    |
    |                          |   0., 0.],                       |
    |                          |         [0., 0., 0., ..., 0.,    |
    |                          |   0., 0.],                       |
    |                          |       ... ...,                   |
    |                          |         [0., 0., 0., ..., 0.,    |
    |                          |   0., 0.],                       |
    |                          |         [0., 0., 0., ..., 0.,    |
    |                          |   0., 0.],                       |
    |                          |         [0., 0., 0., ..., 0.,    |
    |                          |   0., 0.]]])                     |
    | canonical_affine         | [[ 1.19975281e+00                |
    |                          |   -1.89821944e-02                |
    |                          |   -2.72075552e-03                |
    |                          |   -8.58393958e+01]               |
    |                          |  [ 2.45456006e-02                |
    |                          |   9.27821696e-01  1.32986516e-01 |
    |                          |   -1.33765088e+02]               |
    |                          |  [ 5.71511449e-11                |
    |                          |   -1.33014351e-01                |
    |                          |   9.28015888e-01                 |
    |                          |   -8.09051243e+01]               |
    |                          |  [ 0.00000000e+00                |
    |                          |   0.00000000e+00  0.00000000e+00 |
    |                          |    1.00000000e+00]]              |
    +--------------------------+----------------------------------+
    +-----------------------------------+-----------------+
    | Expression                        | Result          |
    +-----------------------------------+-----------------+
    | brain_mri_canonical_data.shape    | (150, 256, 256) |
    | nib.aff2axcodes(canonical_affine) | ('R', 'A', 'S') |
    +-----------------------------------+-----------------+
In [37]:
fig = plt.figure(
    figsize=(figure_size[0], figure_size[1] / 7 * 16), constrained_layout=True
)

gs = gridspec.GridSpec(
    nrows=3, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[49, 50, 50]
)

ax_1 = plt.subplot(gs[0, 0])
ax_2 = plt.subplot(gs[1, 0])
ax_3 = plt.subplot(gs[2, 0])
axs = [ax_1, ax_2, ax_3]

for i in range(3):
    title = (
        f"Canonical NIfTI image grid on the {plane_dict[nib.aff2axcodes(canonical_affine)[i]]} axis "
        f"(the {nth[i+1]} axis)"
    )
    slice_interval, row_size = (
        (8, 8) if brain_mri_canonical_data.shape[i] == 256 else (5, 6)
    )
    if i == 0:
        brain_mri_canonical_swapped = np.swapaxes(brain_mri_canonical_data, 1, 2)
        new_brain_mri_canonical = np.flip(brain_mri_canonical_swapped, 1)
    elif i == 1:
        brain_mri_canonical_swapped = np.swapaxes(brain_mri_canonical_data, 0, 2)
        new_brain_mri_canonical = np.flip(brain_mri_canonical_swapped, 0)
    else:
        brain_mri_canonical_swapped = np.swapaxes(
            brain_mri_canonical_data, axis1=0, axis2=1
        )
        new_brain_mri_canonical = np.flip(brain_mri_canonical_swapped, 1)
    grid_NIfTI_2D_image(
        new_brain_mri_canonical,
        axs[i],
        slice_axis=i,
        slice_interval=slice_interval,
        row_size=row_size,
        title=title,
    )


fig.suptitle(
    "Visual Comparison of Canonical NIfTI Image Grids on Three Different Axes",
    fontsize="x-large",
    x=0.5,
    y=-0.01,
)

plt.show()

Resampling¶

In [38]:
voxel_size = (2, 2, 2)
# `nibabel.processing.conform` samples the given image into a custom shape with custom-sized
# voxels

# The custom size of custom shapes and voxels is determined by the parameters `out_shape` and
# `voxel_size`, which have default values of (256, 256, 256) and 1, where 1 represents the
# size of the voxel in millimeters

# The optional parameter `orientation` specifies the orientation of the output image, and
# its default value is 'RAS'
brain_mri_resized = nibabel.processing.conform(
    brain_mri, out_shape=(128, 128, 100), voxel_size=voxel_size, orientation="PSR"
)
brain_mri_resized_data = brain_mri_resized.get_fdata()
resized_affine = brain_mri_resized.affine

tabulation = Form_Generator()
tabulation.heading_printer("Resampling of NIfTI data")

statements = [
    """
voxel_size = (2, 2, 2)
brain_mri_resized = nibabel.processing.conform(
    brain_mri, out_shape=(128, 128, 100), voxel_size=voxel_size, orientation="PSR"
)
brain_mri_resized_data = brain_mri_resized.get_fdata()
resized_affine = brain_mri_resized.affine
"""
]
tabulation.statement_generator(statements)

variables = ["brain_mri_resized", "brain_mri_resized_data", "resized_affine"]
values = [
    str(brain_mri_resized),
    str(reprlib_rules.repr(brain_mri_resized_data)),
    str(resized_affine),
]
tabulation.variable_generator(variables, values, 3)

expressions = [
    "brain_mri_resized.shape",
    "brain_mri_resized.header.get_zooms()",
    "nib.aff2axcodes(resized_affine)",
    "brain_mri_resized_data.shape",
]

results = [
    str(brain_mri_resized.shape),
    str(brain_mri_resized.header.get_zooms()),
    str(nib.aff2axcodes(resized_affine)),
    str(brain_mri_resized_data.shape),
]
tabulation.expression_generator(expressions, results)
Resampling of NIfTI data

    +--------------------------------------------------------+
    | Statement                                              |
    +--------------------------------------------------------+
    | voxel_size = (2, 2, 2)                                 |
    | brain_mri_resized = nibabel.processing.conform(        |
    |     brain_mri, out_shape=(128, 128, 100),              |
    |     voxel_size=voxel_size, orientation="PSR"           |
    | )                                                      |
    | brain_mri_resized_data = brain_mri_resized.get_fdata() |
    | resized_affine = brain_mri_resized.affine              |
    +--------------------------------------------------------+
    +------------------------+------------------------------------+
    | Variable               | Value                              |
    +------------------------+------------------------------------+
    | brain_mri_resized      | <class                             |
    |                        |    'nibabel.nifti1.Nifti1Image'>   |
    |                        | data shape (128, 128, 100)         |
    |                        | affine:                            |
    |                        | [[ 4.04953482e-02 -5.80427828e-03  |
    |                        |    1.99958157e+00 -9.99979790e+01] |
    |                        |  [-1.97935296e+00  2.83704558e-01  |
    |                        |    4.09092025e-02  1.08523050e+02] |
    |                        |  [ 2.83763951e-01  1.97976717e+00  |
    |                        |    9.52516011e-11 -1.22675404e+02] |
    |                        |  [ 0.00000000e+00  0.00000000e+00  |
    |                        |    0.00000000e+00                  |
    |                        |    1.00000000e+00]]                |
    |                        | metadata:                          |
    |                        | <class                             |
    |                        |    'nibabel.nifti1.Nifti1Header'>  |
    |                        |    object, endian='<'              |
    |                        | sizeof_hdr      : 348              |
    |                        | data_type       : b''              |
    |                        | db_name         : b''              |
    |                        | extents         : 0                |
    |                        | session_error   : 0                |
    |                        | regular         : b'r'             |
    |                        | dim_info        : 0                |
    |                        | dim             : [  3 128 128 100 |
    |                        |      1   1   1   1]                |
    |                        | intent_p1       : 0.0              |
    |                        | intent_p2       : 0.0              |
    |                        | intent_p3       : 0.0              |
    |                        | intent_code     : none             |
    |                        | datatype        : int16            |
    |                        | bitpix          : 16               |
    |                        | slice_start     : 0                |
    |                        | pixdim          : [-1.  2.  2.  2. |
    |                        |     1.  1.  1.  1.]                |
    |                        | vox_offset      : 0.0              |
    |                        | scl_slope       : nan              |
    |                        | scl_inter       : nan              |
    |                        | slice_end       : 0                |
    |                        | slice_code      : unknown          |
    |                        | xyzt_units      : 10               |
    |                        | cal_max         : 0.0              |
    |                        | cal_min         : 0.0              |
    |                        | slice_duration  : 0.0              |
    |                        | toffset         : 0.0              |
    |                        | glmax           : 0                |
    |                        | glmin           : 0                |
    |                        | descrip         : b'MR'            |
    |                        | aux_file        : b''              |
    |                        | qform_code      : unknown          |
    |                        | sform_code      : aligned          |
    |                        | quatern_b       : 0.46861374       |
    |                        | quatern_c       : -0.52952915      |
    |                        | quatern_d       : -0.4576844       |
    |                        | qoffset_x       : -99.99798        |
    |                        | qoffset_y       : 108.52305        |
    |                        | qoffset_z       : -122.67541       |
    |                        | srow_x          : [ 4.0495347e-02  |
    |                        |    -5.8042784e-03  1.9995816e+00   |
    |                        |    -9.9997978e+01]                 |
    |                        | srow_y          : [-1.9793530e+00  |
    |                        |    2.8370455e-01  4.0909201e-02    |
    |                        |    1.0852305e+02]                  |
    |                        | srow_z          : [ 2.8376395e-01  |
    |                        |    1.9797672e+00  9.5251598e-11    |
    |                        |    -1.2267541e+02]                 |
    |                        | intent_name     : b''              |
    |                        | magic           : b'n+1'           |
    | brain_mri_resized_data | array([[[0., 0., 0., ..., 0., 0.,  |
    |                        |    0.],                            |
    |                        |         [0., 0., 0., ..., 0., 0.,  |
    |                        |    0.],                            |
    |                        |         [0., 0., 0., ..., 0., 0.,  |
    |                        |    0.],                            |
    |                        |       ... ...,                     |
    |                        |         [0., 0., 0., ..., 0., 0.,  |
    |                        |    0.],                            |
    |                        |         [0., 0., 0., ..., 0., 0.,  |
    |                        |    0.],                            |
    |                        |         [0., 0., 0., ..., 0., 0.,  |
    |                        |    0.]]])                          |
    | resized_affine         | [[ 4.04953482e-02 -5.80427828e-03  |
    |                        |    1.99958157e+00 -9.99979790e+01] |
    |                        |  [-1.97935296e+00  2.83704558e-01  |
    |                        |    4.09092025e-02  1.08523050e+02] |
    |                        |  [ 2.83763951e-01  1.97976717e+00  |
    |                        |    9.52516011e-11 -1.22675404e+02] |
    |                        |  [ 0.00000000e+00  0.00000000e+00  |
    |                        |    0.00000000e+00                  |
    |                        |    1.00000000e+00]]                |
    +------------------------+------------------------------------+
    +--------------------------------------+-----------------+
    | Expression                           | Result          |
    +--------------------------------------+-----------------+
    | brain_mri_resized.shape              | (128, 128, 100) |
    | brain_mri_resized.header.get_zooms() | (2.0, 2.0, 2.0) |
    | nib.aff2axcodes(resized_affine)      | ('P', 'S', 'R') |
    | brain_mri_resized_data.shape         | (128, 128, 100) |
    +--------------------------------------+-----------------+
In [39]:
fig = plt.figure(
    figsize=(figure_size[0], figure_size[1] / 7 * 16), constrained_layout=True
)

gs = gridspec.GridSpec(
    nrows=3, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[4, 4, 5]
)

ax_1 = plt.subplot(gs[0, 0])
ax_2 = plt.subplot(gs[1, 0])
ax_3 = plt.subplot(gs[2, 0])
axs = [ax_1, ax_2, ax_3]

for i in range(3):
    title = (
        f"Resampled NIfTI image grid on the {plane_dict[nib.aff2axcodes(resized_affine)[i]]} axis "
        f"(the {nth[i+1]} axis)"
    )
    slice_interval, row_size = (
        (4, 8) if brain_mri_resized_data.shape[i] == 128 else (5, 5)
    )
    if i == 0:
        new_brain_mri_resized = np.flip(brain_mri_resized_data, axis=1)
        new_brain_mri_resized = np.flip(new_brain_mri_resized, axis=2)
    elif i == 1:
        new_brain_mri_resized = np.flip(brain_mri_resized_data, 0)
        new_brain_mri_resized = np.flip(new_brain_mri_resized, 2)
    else:
        brain_mri_resized_swapped = np.swapaxes(
            brain_mri_resized_data, axis1=0, axis2=1
        )
        new_brain_mri_resized = np.flip(brain_mri_resized_swapped, 0)
    grid_NIfTI_2D_image(
        new_brain_mri_resized,
        axs[i],
        slice_axis=i,
        slice_interval=slice_interval,
        row_size=row_size,
        title=title,
    )


fig.suptitle(
    "Comparison of the Visualization of the Resampled NIfTI Image Grids on Three Different Axes",
    fontsize="x-large",
    x=0.5,
    y=-0.01,
)

plt.show()
In [ ]: